-
Notifications
You must be signed in to change notification settings - Fork 193
/
incompressible_model.jl
165 lines (141 loc) · 8.41 KB
/
incompressible_model.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
using CUDA: has_cuda
using OrderedCollections: OrderedDict
using Oceananigans: AbstractModel, AbstractOutputWriter, AbstractDiagnostic
using Oceananigans.Architectures: AbstractArchitecture
using Oceananigans.Advection: CenteredSecondOrder
using Oceananigans.Buoyancy: validate_buoyancy, SeawaterBuoyancy
using Oceananigans.BoundaryConditions: regularize_field_boundary_conditions
using Oceananigans.Fields: BackgroundFields, Field, tracernames, VelocityFields, TracerFields, PressureFields
using Oceananigans.Forcings: model_forcing
using Oceananigans.Grids: with_halo
using Oceananigans.Solvers: PressureSolver
using Oceananigans.TimeSteppers: Clock, TimeStepper
using Oceananigans.TurbulenceClosures: ν₀, κ₀, with_tracers, DiffusivityFields, IsotropicDiffusivity
using Oceananigans.Utils: inflate_halo_size, tupleit
mutable struct IncompressibleModel{TS, E, A<:AbstractArchitecture, G, T, B, R, SW, U, C, Φ, F,
V, S, K, BG} <: AbstractModel{TS}
architecture :: A # Computer `Architecture` on which `Model` is run
grid :: G # Grid of physical points on which `Model` is solved
clock :: Clock{T} # Tracks iteration number and simulation time of `Model`
advection :: V # Advection scheme for velocities _and_ tracers
buoyancy :: B # Set of parameters for buoyancy model
coriolis :: R # Set of parameters for the background rotation rate of `Model`
surface_waves :: SW # Set of parameters for surfaces waves via the Craik-Leibovich approximation
forcing :: F # Container for forcing functions defined by the user
closure :: E # Diffusive 'turbulence closure' for all model fields
background_fields :: BG # Background velocity and tracer fields
velocities :: U # Container for velocity fields `u`, `v`, and `w`
tracers :: C # Container for tracer fields
pressures :: Φ # Container for hydrostatic and nonhydrostatic pressure
diffusivities :: K # Container for turbulent diffusivities
timestepper :: TS # Object containing timestepper fields and parameters
pressure_solver :: S # Pressure/Poisson solver
end
"""
IncompressibleModel(;
grid,
architecture = CPU(),
float_type = Float64,
clock = Clock{float_type}(0, 0, 1),
advection = CenteredSecondOrder(),
buoyancy = SeawaterBuoyancy(float_type),
coriolis = nothing,
surface_waves = nothing,
forcing = NamedTuple(),
closure = IsotropicDiffusivity(float_type, ν=ν₀, κ=κ₀),
boundary_conditions = NamedTuple(),
tracers = (:T, :S),
timestepper = :QuasiAdamsBashforth2,
background_fields = NamedTuple(),
velocities = nothing,
pressures = nothing,
diffusivities = nothing,
pressure_solver = nothing
)
Construct an incompressible `Oceananigans.jl` model on `grid`.
Keyword arguments
=================
- `grid`: (required) The resolution and discrete geometry on which `model` is solved.
- `architecture`: `CPU()` or `GPU()`. The computer architecture used to time-step `model`.
- `float_type`: `Float32` or `Float64`. The floating point type used for `model` data.
- `advection`: The scheme that advects velocities and tracers. See `Oceananigans.Advection`.
- `buoyancy`: The buoyancy model. See `Oceananigans.Buoyancy`.
- `closure`: The turbulence closure for `model`. See `Oceananigans.TurbulenceClosures`.
- `coriolis`: Parameters for the background rotation rate of the model.
- `forcing`: `NamedTuple` of user-defined forcing functions that contribute to solution tendencies.
- `boundary_conditions`: `NamedTuple` containing field boundary conditions.
- `tracers`: A tuple of symbols defining the names of the modeled tracers, or a `NamedTuple` of
preallocated `CellField`s.
- `timestepper`: A symbol that specifies the time-stepping method. Either `:QuasiAdamsBashforth2` or
`:RungeKutta3`.
"""
function IncompressibleModel(;
grid,
architecture::AbstractArchitecture = CPU(),
float_type = Float64,
clock = Clock{float_type}(0, 0, 1),
advection = CenteredSecondOrder(),
buoyancy = SeawaterBuoyancy(float_type),
coriolis = nothing,
surface_waves = nothing,
forcing::NamedTuple = NamedTuple(),
closure = IsotropicDiffusivity(float_type, ν=ν₀, κ=κ₀),
boundary_conditions::NamedTuple = NamedTuple(),
tracers = (:T, :S),
timestepper = :QuasiAdamsBashforth2,
background_fields::NamedTuple = NamedTuple(),
velocities = nothing,
pressures = nothing,
diffusivities = nothing,
pressure_solver = nothing
)
if architecture == GPU() && !has_cuda()
throw(ArgumentError("Cannot create a GPU model. No CUDA-enabled GPU was detected!"))
end
tracers = tupleit(tracers) # supports tracers=:c keyword argument (for example)
validate_buoyancy(buoyancy, tracernames(tracers))
# Adjust halos when the advection scheme or turbulence closure requires it.
# Note that halos are isotropic by default; however we respect user-input here
# by adjusting each (x, y, z) halo individually.
Hx, Hy, Hz = inflate_halo_size(grid.Hx, grid.Hy, grid.Hz, advection, closure)
grid = with_halo((Hx, Hy, Hz), grid)
# Recursively "regularize" field-dependent boundary conditions by supplying list of tracer names.
# We also regularize boundary conditions included in velocities, tracers, pressures, and diffusivities.
# Note that we do not regularize boundary conditions contained in *tupled* diffusivity fields right now.
embedded_boundary_conditions = merge(extract_boundary_conditions(velocities),
extract_boundary_conditions(tracers),
extract_boundary_conditions(pressures),
extract_boundary_conditions(diffusivities))
boundary_conditions = merge(embedded_boundary_conditions, boundary_conditions)
boundary_conditions = regularize_field_boundary_conditions(boundary_conditions, grid, tracernames(tracers), nothing)
# Either check grid-correctness, or construct tuples of fields
velocities = VelocityFields(velocities, architecture, grid, boundary_conditions)
tracers = TracerFields(tracers, architecture, grid, boundary_conditions)
pressures = PressureFields(pressures, architecture, grid, boundary_conditions)
diffusivities = DiffusivityFields(diffusivities, architecture, grid,
tracernames(tracers), boundary_conditions, closure)
pressure_solver = PressureSolver(pressure_solver, architecture, grid, PressureBoundaryConditions(grid))
background_fields = BackgroundFields(background_fields, tracernames(tracers), grid, clock)
# Instantiate timestepper if not already instantiated
timestepper = TimeStepper(timestepper, architecture, grid, tracernames(tracers))
# Regularize forcing and closure for model tracer and velocity fields.
forcing = model_forcing(tracernames(tracers); forcing...)
closure = with_tracers(tracernames(tracers), closure)
return IncompressibleModel(architecture, grid, clock, advection, buoyancy, coriolis, surface_waves,
forcing, closure, background_fields, velocities, tracers, pressures,
diffusivities, timestepper, pressure_solver)
end
#####
##### Recursive util for building NamedTuples of boundary conditions from NamedTuples of fields
#####
##### Note: ignores tuples, including tuples of Symbols (tracer names) and
##### tuples of DiffusivityFields (which occur for tupled closures)
#####
extract_boundary_conditions(::Nothing) = NamedTuple()
extract_boundary_conditions(::Tuple) = NamedTuple()
function extract_boundary_conditions(field_tuple::NamedTuple)
names = propertynames(field_tuple)
bcs = Tuple(extract_boundary_conditions(field) for field in field_tuple)
return NamedTuple{names}(bcs)
end
extract_boundary_conditions(field::Field) = field.boundary_conditions