diff --git a/.github/workflows/Whitespace.yml b/.github/workflows/Whitespace.yml new file mode 100644 index 00000000000..7059b8775da --- /dev/null +++ b/.github/workflows/Whitespace.yml @@ -0,0 +1,36 @@ +name: Whitespace + +permissions: {} + +on: + push: + branches: + - 'main' + pull_request: + +jobs: + whitespace: + name: Check whitespace + runs-on: ubuntu-latest + timeout-minutes: 2 + steps: + - name: Checkout the repository + uses: actions/checkout@v6 + with: + persist-credentials: false + - name: Checkout the JuliaLang/julia repository + uses: actions/checkout@v6 + with: + persist-credentials: false + repository: 'JuliaLang/julia' + # Clone Julia in a subdir. + path: '.julia' + # Check out a fixed revision to avoid surprises in case the script is + # changed in the future. + ref: '3b12a882e887753e4d2e9e9db65d99d3f7d9e26b' + - uses: julia-actions/setup-julia@v2 + with: + version: '1.12.4' + - name: Check whitespace + run: | + .julia/contrib/check-whitespace.jl diff --git a/AGENTS.md b/AGENTS.md index 11db4d83c55..d3c99ce29a6 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -5,7 +5,7 @@ Oceananigans.jl is a Julia package for fast, friendly, flexible, ocean-flavored fluid dynamics simulations on CPUs and GPUs. It provides a framework for solving the incompressible (or Boussinesq) Navier-Stokes equations with various model configurations including: - Nonhydrostatic models with free surfaces -- Hydrostatic models for large-scale ocean simulations +- Hydrostatic models for large-scale ocean simulations - Shallow water models - Support for a variety of grids: RectilinearGrid, LatitudeLongitudeGrid, CubedSphereGrid - Support for complex domains using ImmersedBoundaryGrid @@ -22,10 +22,10 @@ It provides a framework for solving the incompressible (or Boussinesq) Navier-St 1. **Explicit Imports**: Use `ExplicitImports.jl` style - explicitly import all used functions/types - Import from modules explicitly (already done in src/Oceananigans.jl) - Tests automatically check for proper imports - + 2. **Type Stability**: Prioritize type-stable code for performance - All structs must be concretely typed - + 3. **Kernel Functions**: For GPU compatibility: - Use KernelAbstractions.jl syntax for kernels, eg `@kernel`, `@index` - Keep kernels type-stable and allocation-free @@ -35,7 +35,7 @@ It provides a framework for solving the incompressible (or Boussinesq) Navier-St - Mark functions called inside kernels with `@inline` - **Never use loops outside kernels**: Always replace `for` loops that iterate over grid points with kernels launched via `launch!`. This ensures code works on both CPU and GPU. - + 4. **Documentation**: - Use DocStringExtensions.jl for consistent docstrings - Include `$(SIGNATURES)` for automatic signature documentation @@ -53,7 +53,7 @@ It provides a framework for solving the incompressible (or Boussinesq) Navier-St - **NonhydrostaticModel**: `NonhydrostaticModel(grid; ...)` - `grid` is positional - **ShallowWaterModel**: `ShallowWaterModel(grid; gravitational_acceleration, ...)` - both `grid` and `gravitational_acceleration` are positional - **Important**: When there are no keyword arguments, omit the semicolon: - - ✅ `NonhydrostaticModel(grid)` + - ✅ `NonhydrostaticModel(grid)` - ❌ `NonhydrostaticModel(grid;)` - When keyword arguments are present, use the semicolon: - ✅ `NonhydrostaticModel(grid; closure=nothing)` @@ -350,7 +350,7 @@ serve(dir="docs/build") When implementing a simulation from a published paper: ### 1. Parameter Extraction -- **Read the paper carefully** and extract ALL parameters: domain size, resolution, physical constants, +- **Read the paper carefully** and extract ALL parameters: domain size, resolution, physical constants, boundary conditions, initial conditions, forcing, closure parameters - Look for parameter tables (often "Table 1" or similar) - Check figure captions for additional details @@ -393,11 +393,11 @@ Before running a long simulation: - Quantitative comparison: compute the same diagnostics as the paper ### 7. Common Issues -- **NaN blowups**: Usually from timestep too large, unstable initial conditions, +- **NaN blowups**: Usually from timestep too large, unstable initial conditions, or if-else statements on GPU (use `ifelse` instead) -- **Nothing happening**: Check that buoyancy anomaly has the right sign, +- **Nothing happening**: Check that buoyancy anomaly has the right sign, that initial conditions are actually applied, that forcing is active -- **Wrong direction of flow**: Check coordinate conventions (is y increasing +- **Wrong direction of flow**: Check coordinate conventions (is y increasing upslope or downslope?) - **GPU issues**: Avoid branching, ensure type stability, use `randn()` carefully diff --git a/benchmark/README.md b/benchmark/README.md index 1e007053005..4f376f9680f 100644 --- a/benchmark/README.md +++ b/benchmark/README.md @@ -23,4 +23,3 @@ Run distributed benchmarks by running the launcher scripts for either the shallo ## Measuring performance regression Running the `benchmark_regression.jl` script will run the nonhydrostatic model tests on the current branch and on the main branch for comparison. This is useful to test whether the current branch slows down the code or introduces any performance regression. - diff --git a/benchmark/benchmark_regression.jl b/benchmark/benchmark_regression.jl index e99502b39dd..1da7cbb9af6 100644 --- a/benchmark/benchmark_regression.jl +++ b/benchmark/benchmark_regression.jl @@ -17,4 +17,3 @@ for (case, trial) in results println("Results for $case") display(trial) end - diff --git a/benchmark/benchmarkable_nonhydrostatic_model.jl b/benchmark/benchmarkable_nonhydrostatic_model.jl index f7284a4b71b..e1a19861431 100644 --- a/benchmark/benchmarkable_nonhydrostatic_model.jl +++ b/benchmark/benchmarkable_nonhydrostatic_model.jl @@ -28,4 +28,3 @@ for Arch in Architectures, N in Ns SUITE[(Arch, N)] = benchmark end - diff --git a/docs/src/developer_docs/model_interface.md b/docs/src/developer_docs/model_interface.md index 1c1cb47f6e2..25d6c867e4c 100644 --- a/docs/src/developer_docs/model_interface.md +++ b/docs/src/developer_docs/model_interface.md @@ -65,17 +65,17 @@ implement (or inherit sane fallbacks for) the items listed below. - `update_state!(model, callbacks=[]; compute_tendencies=true)`: invoked by Simulation right after `initialize!` and inside most time steppers. This is where models fill halos, update boundary conditions, recompute auxiliary - fields, and run [`Callback`](@ref callbacks)s with an `UpdateStateCallsite`. - PDE-based models typically finish by calling `compute_tendencies!(model, callbacks)` + fields, and run [`Callback`](@ref callbacks)s with an `UpdateStateCallsite`. + PDE-based models typically finish by calling `compute_tendencies!(model, callbacks)` so that any `TendencyCallsite` callbacks can modify tendencies before integration. Note that `compute_tendencies!` is not part of the required interface—it is simply a useful pattern for models that integrate differential equations. - `time_step!(model, Δt; callbacks=[])`: advances the model clock and its prognostic variables by one step. Simulation hands in the tuple of - `ModelCallsite` [`Callback`](@ref callbacks)s so the model can execute - `TendencyCallsite` (before tendencies are applied) and `UpdateStateCallsite` - callbacks (after auxiliary updates). The method must call `tick!(model.clock, Δt)` + `ModelCallsite` [`Callback`](@ref callbacks)s so the model can execute + `TendencyCallsite` (before tendencies are applied) and `UpdateStateCallsite` + callbacks (after auxiliary updates). The method must call `tick!(model.clock, Δt)` (or equivalent) so that `time(model)` and `iteration(model)` remain consistent. - `set!(model, kw...)`: not strictly required, but strongly recommended as an @@ -241,8 +241,8 @@ fig This minimal implementation inherits all other behavior from the generic `AbstractModel` fallbacks: Simulation can query `time(sim.model)`, diagnostics -can read `sim.model.clock`, and [`Callback`](@ref callbacks)s scheduled on -`ModelCallsite`s execute because `time_step!` forwards the tuple that Simulation +can read `sim.model.clock`, and [`Callback`](@ref callbacks)s scheduled on +`ModelCallsite`s execute because `time_step!` forwards the tuple that Simulation hands to it. Note that this model has no grid, no fields, and no time-stepper object—just the essentials. Larger models can follow the same recipe while adding grids, fields, closures, and time steppers as needed. @@ -411,4 +411,3 @@ operators within a custom `AbstractModel`. The key additions compared to the - Third-order Runge-Kutta (RK3) time-stepping using Williamson's low-storage scheme - Using `fill_halo_regions!` in `update_state!` to maintain periodic boundary conditions - Leveraging `AbstractOperations` (`∂x`) for computing spatial derivatives via broadcasting - diff --git a/docs/src/developer_docs/turbulence_closures.md b/docs/src/developer_docs/turbulence_closures.md index 2962ffc3cd2..5f07eee6462 100644 --- a/docs/src/developer_docs/turbulence_closures.md +++ b/docs/src/developer_docs/turbulence_closures.md @@ -234,7 +234,7 @@ Good display methods help users understand their closures: ```@example pp_closure using Oceananigans.Utils: prettysummary -Base.summary(closure::PPVD{TD}) where TD = +Base.summary(closure::PPVD{TD}) where TD = string("PacanowskiPhilanderVerticalDiffusivity{$TD}") function Base.show(io::IO, closure::PPVD) @@ -277,7 +277,7 @@ We'll create a helper function to set up and run simulations: ```@example pp_closure function run_boundary_layer(closure; stop_time) grid = RectilinearGrid(size=Nz, z=(-Lz, 0), topology=(Flat, Flat, Bounded)) - + u_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(τˣ)) b_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Jᵇ)) @@ -287,13 +287,13 @@ function run_boundary_layer(closure; stop_time) tracers = :b, coriolis = FPlane(f=f), boundary_conditions = (u=u_bcs, b=b_bcs)) - + set!(model, b = z -> N² * z) # linear stratification - + simulation = Simulation(model; Δt=10minutes, stop_time) conjure_time_step_wizard!(simulation, cfl=0.5, max_Δt=10minutes) run!(simulation) - + return model end nothing # hide @@ -331,22 +331,22 @@ z_tked = znodes(model_tked.tracers.b) ax1 = Axis(fig[1, 1], xlabel="Buoyancy (m s⁻²)", ylabel="z (m)", title="Buoyancy profile") -lpp = lines!(ax1, interior(model_pp.tracers.b, 1, 1, :), z_pp, +lpp = lines!(ax1, interior(model_pp.tracers.b, 1, 1, :), z_pp, label="Pacanowski-Philander", linewidth=2) -lcatke = lines!(ax1, interior(model_catke.tracers.b, 1, 1, :), z_catke, +lcatke = lines!(ax1, interior(model_catke.tracers.b, 1, 1, :), z_catke, label="CATKE", linewidth=2, linestyle=:dash) -ltked = lines!(ax1, interior(model_tked.tracers.b, 1, 1, :), z_tked, +ltked = lines!(ax1, interior(model_tked.tracers.b, 1, 1, :), z_tked, label="TKE-Dissipation", linewidth=2, linestyle=:dot) ## Velocity profiles ax2 = Axis(fig[1, 2], xlabel="Velocity (m s⁻¹)", ylabel="z (m)", title="Zonal velocity") -lines!(ax2, interior(model_pp.velocities.u, 1, 1, :), z_pp, +lines!(ax2, interior(model_pp.velocities.u, 1, 1, :), z_pp, label="PP", linewidth=2) -lines!(ax2, interior(model_catke.velocities.u, 1, 1, :), z_catke, +lines!(ax2, interior(model_catke.velocities.u, 1, 1, :), z_catke, label="CATKE", linewidth=2, linestyle=:dash) -lines!(ax2, interior(model_tked.velocities.u, 1, 1, :), z_tked, +lines!(ax2, interior(model_tked.velocities.u, 1, 1, :), z_tked, label="TKE-ϵ", linewidth=2, linestyle=:dot) ## Diffusivity profiles @@ -472,7 +472,7 @@ If you'd like to contribute your closure to Oceananigans itself, here are the ad ### 1. Create a source file -Place your implementation in a file under +Place your implementation in a file under `src/TurbulenceClosures/turbulence_closure_implementations/`. For example: ``` diff --git a/docs/src/field_time_series.md b/docs/src/field_time_series.md index a9dfb0cf8a8..49811ad2ec1 100644 --- a/docs/src/field_time_series.md +++ b/docs/src/field_time_series.md @@ -1,2 +1 @@ # FieldTimeSeries - diff --git a/docs/src/gallery.md b/docs/src/gallery.md index 71c319b3967..af960e6d5ce 100644 --- a/docs/src/gallery.md +++ b/docs/src/gallery.md @@ -63,4 +63,3 @@ the fine details at the surface as it cools and the layers of different temperat by internal waves. [![Watch free convection with wind stress in action](https://raw.githubusercontent.com/ali-ramadhan/ali-ramadhan.Github.io/master/img/wind_stress_unstable_7500.png)](https://www.youtube.com/watch?v=ob6OMQgPfI4) - diff --git a/docs/src/models/stokes_drift.md b/docs/src/models/stokes_drift.md index 9c0c4794ad7..e803fd41429 100644 --- a/docs/src/models/stokes_drift.md +++ b/docs/src/models/stokes_drift.md @@ -207,4 +207,3 @@ The function signature for [`StokesDrift`](@ref) depends on the grid topology: functions are called as `f(x, z, t)` (the ``y`` coordinate is omitted). - When `parameters` is provided, it is passed as an additional final argument: `f(x, y, z, t, parameters)` or `f(x, z, t, parameters)`. - diff --git a/docs/src/numerical_implementation/elliptic_solvers.md b/docs/src/numerical_implementation/elliptic_solvers.md index 66ad307550a..cfdb0b6c32e 100644 --- a/docs/src/numerical_implementation/elliptic_solvers.md +++ b/docs/src/numerical_implementation/elliptic_solvers.md @@ -16,7 +16,7 @@ along with homogenous Neumann boundary conditions ``\boldsymbol{v} \cdot \boldsy denotes the source term for the Poisson equation. !!! note "Dividing small timesteps" - In practice, in order to avoid division by extremely small numbers when ``\Delta t \lesssim \epsilon``, + In practice, in order to avoid division by extremely small numbers when ``\Delta t \lesssim \epsilon``, we solve the Poisson equation for ``p_{NH} \Delta t`` instead. ``\Delta t`` is then removed from the pressure field after the elliptic solve routine. diff --git a/docs/src/operations.md b/docs/src/operations.md index 4fa6fdbad95..5baae4a10b3 100644 --- a/docs/src/operations.md +++ b/docs/src/operations.md @@ -43,8 +43,8 @@ BinaryOperation at (Center, Center, Center) ├── grid: 4×1×4 RectilinearGrid{Float64, Periodic, Flat, Bounded} on CPU with 3×0×3 halo └── tree: * at (Center, Center, Center) -    ├── 2 -    └── 4×1×4 Field{Center, Center, Center} on RectilinearGrid on CPU + ├── 2 + └── 4×1×4 Field{Center, Center, Center} on RectilinearGrid on CPU ``` and even by chaining expressions together, which may themselves include `AbstractOperations`, @@ -57,13 +57,13 @@ MultiaryOperation at (Center, Center, Center) ├── grid: 4×1×4 RectilinearGrid{Float64, Periodic, Flat, Bounded} on CPU with 3×0×3 halo └── tree: + at (Center, Center, Center) -    ├── ^ at (Center, Center, Center) -    │   ├── 4×1×4 Field{Center, Center, Center} on RectilinearGrid on CPU -    │   └── 2 -    ├── * at (Center, Center, Center) -    │   ├── 2 -    │   └── 4×1×4 Field{Center, Center, Center} on RectilinearGrid on CPU -    └── 1 + ├── ^ at (Center, Center, Center) + │ ├── 4×1×4 Field{Center, Center, Center} on RectilinearGrid on CPU + │ └── 2 + ├── * at (Center, Center, Center) + │ ├── 2 + │ └── 4×1×4 Field{Center, Center, Center} on RectilinearGrid on CPU + └── 1 ``` Like `Field`s, `AbstractOperations` have a location and a grid. @@ -77,7 +77,7 @@ UnaryOperation at (Center, Center, Center) ├── grid: 4×1×4 RectilinearGrid{Float64, Periodic, Flat, Bounded} on CPU with 3×0×3 halo └── tree: cos at (Center, Center, Center) via identity -    └── 4×1×4 Field{Center, Center, Center} on RectilinearGrid on CPU + └── 4×1×4 Field{Center, Center, Center} on RectilinearGrid on CPU ``` ```jldoctest operations @@ -88,10 +88,10 @@ MultiaryOperation at (Center, Center, Center) ├── grid: 4×1×4 RectilinearGrid{Float64, Periodic, Flat, Bounded} on CPU with 3×0×3 halo └── tree: + at (Center, Center, Center) -    ├── 4×1×4 Field{Center, Center, Center} on RectilinearGrid on CPU -    ├── 4×1×4 Field{Center, Center, Center} on RectilinearGrid on CPU -    ├── 4×1×4 Field{Center, Center, Center} on RectilinearGrid on CPU -    └── 4×1×4 Field{Center, Center, Center} on RectilinearGrid on CPU + ├── 4×1×4 Field{Center, Center, Center} on RectilinearGrid on CPU + ├── 4×1×4 Field{Center, Center, Center} on RectilinearGrid on CPU + ├── 4×1×4 Field{Center, Center, Center} on RectilinearGrid on CPU + └── 4×1×4 Field{Center, Center, Center} on RectilinearGrid on CPU ``` `UnaryOperation`, `BinaryOperation` and `MultiaryOperation` all have both an "operator", and between 1 and many. @@ -105,7 +105,7 @@ Derivative at (Face, Center, Center) ├── grid: 4×1×4 RectilinearGrid{Float64, Periodic, Flat, Bounded} on CPU with 3×0×3 halo └── tree: ∂xᶠᶜᶜ at (Face, Center, Center) via identity -    └── 4×1×4 Field{Center, Center, Center} on RectilinearGrid on CPU + └── 4×1×4 Field{Center, Center, Center} on RectilinearGrid on CPU ``` !!! note @@ -158,14 +158,14 @@ BinaryOperation at (Face, Center, Center) ├── grid: 4×1×4 RectilinearGrid{Float64, Periodic, Flat, Bounded} on CPU with 3×0×3 halo └── tree: + at (Face, Center, Center) -    ├── ^ at (Face, Center, Center) -    │   ├── ∂xᶠᶜᶜ at (Face, Center, Center) via identity -    │   │   └── 4×1×4 Field{Center, Center, Center} on RectilinearGrid on CPU -    │   └── 2 -    └── ^ at (Center, Center, Face) -       ├── ∂zᶜᶜᶠ at (Center, Center, Face) via identity -       │   └── 4×1×4 Field{Center, Center, Center} on RectilinearGrid on CPU -       └── 2 + ├── ^ at (Face, Center, Center) + │ ├── ∂xᶠᶜᶜ at (Face, Center, Center) via identity + │ │ └── 4×1×4 Field{Center, Center, Center} on RectilinearGrid on CPU + │ └── 2 + └── ^ at (Center, Center, Face) + ├── ∂zᶜᶜᶠ at (Center, Center, Face) via identity + │ └── 4×1×4 Field{Center, Center, Center} on RectilinearGrid on CPU + └── 2 ``` Because `∂x(c)^2` is located at `(Face, Center, Center)` and `∂z(c)^2` is located at `(Center, Center, Face)`, @@ -181,14 +181,14 @@ BinaryOperation at (Center, Center, Face) ├── grid: 4×1×4 RectilinearGrid{Float64, Periodic, Flat, Bounded} on CPU with 3×0×3 halo └── tree: + at (Center, Center, Face) -    ├── ^ at (Center, Center, Face) -    │   ├── ∂zᶜᶜᶠ at (Center, Center, Face) via identity -    │   │   └── 4×1×4 Field{Center, Center, Center} on RectilinearGrid on CPU -    │   └── 2 -    └── ^ at (Face, Center, Center) -       ├── ∂xᶠᶜᶜ at (Face, Center, Center) via identity -       │   └── 4×1×4 Field{Center, Center, Center} on RectilinearGrid on CPU -       └── 2 + ├── ^ at (Center, Center, Face) + │ ├── ∂zᶜᶜᶠ at (Center, Center, Face) via identity + │ │ └── 4×1×4 Field{Center, Center, Center} on RectilinearGrid on CPU + │ └── 2 + └── ^ at (Face, Center, Center) + ├── ∂xᶠᶜᶜ at (Face, Center, Center) via identity + │ └── 4×1×4 Field{Center, Center, Center} on RectilinearGrid on CPU + └── 2 ``` ends up at `(Center, Center, Face)`. To control the location of an operation we use the macro `@at`, @@ -201,14 +201,14 @@ BinaryOperation at (Center, Center, Center) ├── grid: 4×1×4 RectilinearGrid{Float64, Periodic, Flat, Bounded} on CPU with 3×0×3 halo └── tree: + at (Center, Center, Center) -    ├── ^ at (Center, Center, Center) -    │   ├── ∂xᶠᶜᶜ at (Center, Center, Center) via ℑxᶜᵃᵃ -    │   │   └── 4×1×4 Field{Center, Center, Center} on RectilinearGrid on CPU -    │   └── 2 -    └── ^ at (Center, Center, Center) -       ├── ∂zᶜᶜᶠ at (Center, Center, Center) via ℑzᵃᵃᶜ -       │   └── 4×1×4 Field{Center, Center, Center} on RectilinearGrid on CPU -       └── 2 + ├── ^ at (Center, Center, Center) + │ ├── ∂xᶠᶜᶜ at (Center, Center, Center) via ℑxᶜᵃᵃ + │ │ └── 4×1×4 Field{Center, Center, Center} on RectilinearGrid on CPU + │ └── 2 + └── ^ at (Center, Center, Center) + ├── ∂zᶜᶜᶠ at (Center, Center, Center) via ℑzᵃᵃᶜ + │ └── 4×1×4 Field{Center, Center, Center} on RectilinearGrid on CPU + └── 2 ``` ## Averages and integrals diff --git a/docs/src/physics/boussinesq.md b/docs/src/physics/boussinesq.md index 2f1b9fa6615..cb3d2db3787 100644 --- a/docs/src/physics/boussinesq.md +++ b/docs/src/physics/boussinesq.md @@ -50,4 +50,3 @@ abuse notation a bit and denote the kinematic pressure simply as ``p``. for an oceanographic introduction to the Boussinesq equations and Vallis (2017, Section 2.A) for an asymptotic derivation. See Kundu (2015, Section 4.9) for an engineering introduction. - diff --git a/ext/OceananigansEnzymeExt.jl b/ext/OceananigansEnzymeExt.jl index faca70c2b8e..ed623abd6c7 100644 --- a/ext/OceananigansEnzymeExt.jl +++ b/ext/OceananigansEnzymeExt.jl @@ -133,8 +133,8 @@ EnzymeCore.EnzymeRules.inactive(::typeof(Oceananigans.BoundaryConditions.periodi # # FunctionFieldType(dfunction_field_func, grid.val; clock, parameters) # else -# ntuple(Val(config_width)) do i -# Base.@_inline_meta +# ntuple(Val(config_width)) do i +# Base.@_inline_meta # # dfunction_field_func = if function_field_is_active # dactives[i][] @@ -143,7 +143,7 @@ EnzymeCore.EnzymeRules.inactive(::typeof(Oceananigans.BoundaryConditions.periodi # end # # FunctionFieldType(dfunction_field_func, grid.val; clock, parameters) -# end +# end # end # # P = EnzymeCore.EnzymeRules.needs_primal(config) ? RT : Nothing diff --git a/ext/OceananigansMetalExt.jl b/ext/OceananigansMetalExt.jl index 6cc0ce177d6..80f03ed9976 100644 --- a/ext/OceananigansMetalExt.jl +++ b/ext/OceananigansMetalExt.jl @@ -51,4 +51,3 @@ Metal.@device_override @inline function __validindex(ctx::MappedCompilerMetadata end end # module - diff --git a/ext/OceananigansReactantExt/Grids/Grids.jl b/ext/OceananigansReactantExt/Grids/Grids.jl index c1b3106cf8f..bb0373a29da 100644 --- a/ext/OceananigansReactantExt/Grids/Grids.jl +++ b/ext/OceananigansReactantExt/Grids/Grids.jl @@ -90,4 +90,3 @@ function offset_data(underlying_data::ConcreteRArray, loc, topo, N, H, indices:: end end # module - diff --git a/ext/OceananigansReactantExt/Grids/sharded_grids.jl b/ext/OceananigansReactantExt/Grids/sharded_grids.jl index 20f2351bb1b..2f67e273517 100644 --- a/ext/OceananigansReactantExt/Grids/sharded_grids.jl +++ b/ext/OceananigansReactantExt/Grids/sharded_grids.jl @@ -263,4 +263,3 @@ function Oceananigans.Grids.zeros(arch::ShardedDistributed, FT, global_sz...) sharding=Sharding.DimsSharding(arch.connectivity, (1, 2), (:x, :y)) ) end - diff --git a/ext/OceananigansReactantExt/OceananigansReactantExt.jl b/ext/OceananigansReactantExt/OceananigansReactantExt.jl index 6e832e06a08..e46592ae5ff 100644 --- a/ext/OceananigansReactantExt/OceananigansReactantExt.jl +++ b/ext/OceananigansReactantExt/OceananigansReactantExt.jl @@ -134,10 +134,10 @@ Base.@nospecializeinfer function Reactant.traced_type_inner( CF2 = Reactant.traced_type_inner(CF, seen, mode, track_numbers, sharding, runtime) FF2 = Reactant.traced_type_inner(FF, seen, mode, track_numbers, sharding, runtime) for NF in (CC2, FC2, CF2, FF2) - if NF === Nothing - continue - end - FT2 = Reactant.promote_traced_type(FT2, eltype(NF)) + if NF === Nothing + continue + end + FT2 = Reactant.promote_traced_type(FT2, eltype(NF)) end rFT2 = Reactant.traced_type_inner(rFT, seen, mode, track_numbers, sharding, runtime) return Oceananigans.Grids.OrthogonalSphericalShellGrid{FT2, TX2, TY2, TZ2, Z2, Map2, CC2, FC2, CF2, FF2, Arch, rFT2} @@ -171,14 +171,14 @@ Base.@nospecializeinfer function Reactant.traced_type_inner( end Base.@nospecializeinfer function Reactant.traced_type_inner( - @nospecialize(OA::Type{LatitudeLongitudeGrid{FT, TX, TY, TZ, Z, DXF, DXC, XF, XC, DYF, DYC, YF, YC, + @nospecialize(OA::Type{LatitudeLongitudeGrid{FT, TX, TY, TZ, Z, DXF, DXC, XF, XC, DYF, DYC, YF, YC, DXCC, DXFC, DXCF, DXFF, DYFC, DYCF, Arch, I}}), seen, mode::Reactant.TraceMode, @nospecialize(track_numbers::Type), @nospecialize(sharding), @nospecialize(runtime) -) where {FT, TX, TY, TZ, Z, DXF, DXC, XF, XC, DYF, DYC, YF, YC, DXCC, DXFC, DXCF, DXFF, DYFC, DYCF, Arch, I} +) where {FT, TX, TY, TZ, Z, DXF, DXC, XF, XC, DYF, DYC, YF, YC, DXCC, DXFC, DXCF, DXFF, DYFC, DYCF, Arch, I} TX2 = Reactant.traced_type_inner(TX, seen, mode, track_numbers, sharding, runtime) TY2 = Reactant.traced_type_inner(TY, seen, mode, track_numbers, sharding, runtime) TZ2 = Reactant.traced_type_inner(TZ, seen, mode, track_numbers, sharding, runtime) @@ -202,14 +202,14 @@ Base.@nospecializeinfer function Reactant.traced_type_inner( FT2 = Reactant.traced_type_inner(FT, seen, mode, track_numbers, sharding, runtime) for NF in (XF2, XC2, YF2, YC2, DXCC2, DXFC2, DYCF2, DYCF2, DXFF2) - if NF === Nothing - continue - end - FT2 = Reactant.promote_traced_type(FT2, eltype(NF)) + if NF === Nothing + continue + end + FT2 = Reactant.promote_traced_type(FT2, eltype(NF)) end - res = Oceananigans.Grids.LatitudeLongitudeGrid{FT2, TX2, TY2, TZ2, Z2, DXF2, DXC2, XF2, XC2, DYF2, DYC2, YF2, YC2, - DXCC2, DXFC2, DXCF2, DXFF2, DYFC2, DYCF2, Arch, I2} + res = Oceananigans.Grids.LatitudeLongitudeGrid{FT2, TX2, TY2, TZ2, Z2, DXF2, DXC2, XF2, XC2, DYF2, DYC2, YF2, YC2, + DXCC2, DXFC2, DXCF2, DXFF2, DYFC2, DYCF2, Arch, I2} return res end @@ -266,13 +266,13 @@ end @assert size(tvals) == size(c) gf = Reactant.call_with_reactant(getindex, c.operand, axes2...) if gf isa AbstractFloat - gf = Reactant.Ops.fill(gf, size(c)) + gf = Reactant.Ops.fill(gf, size(c)) end Reactant.TracedRArrayOverrides._copyto!(tvals, Base.broadcasted(c.func isa Nothing ? Base.identity : c.func, gf)) mask = c.mask if mask isa AbstractFloat && typeof(mask) != Reactant.unwrapped_eltype(Base.eltype(c)) - mask = Base.eltype(c)(mask) + mask = Base.eltype(c)(mask) end return Reactant.Ops.select( @@ -348,4 +348,3 @@ Base.getindex(array::OffsetVector{T, <:Reactant.AbstractConcreteArray{T, 1}}, :: # using .Solvers end # module - diff --git a/ext/OceananigansReactantExt/OutputReaders.jl b/ext/OceananigansReactantExt/OutputReaders.jl index 1259ef29e01..e28e50f8994 100644 --- a/ext/OceananigansReactantExt/OutputReaders.jl +++ b/ext/OceananigansReactantExt/OutputReaders.jl @@ -23,7 +23,7 @@ import Oceananigans.OutputReaders: find_time_index, cpu_interpolating_time_indic return ñ, n₁, n₂ end -function cpu_interpolating_time_indices(::ReactantState, times, time_indexing, t) +function cpu_interpolating_time_indices(::ReactantState, times, time_indexing, t) cpu_times = on_architecture(Oceananigans.Architectures.CPU(), times) return TimeInterpolator(time_indexing, times, t) end diff --git a/ext/OceananigansReactantExt/Simulations/run.jl b/ext/OceananigansReactantExt/Simulations/run.jl index 02e48b44485..21f39fdd1d0 100644 --- a/ext/OceananigansReactantExt/Simulations/run.jl +++ b/ext/OceananigansReactantExt/Simulations/run.jl @@ -12,4 +12,3 @@ function time_step_for!(sim::ReactantSimulation, Nsteps) end return nothing end - diff --git a/ext/OceananigansReactantExt/Simulations/simulation.jl b/ext/OceananigansReactantExt/Simulations/simulation.jl index 67709114113..543823de61e 100644 --- a/ext/OceananigansReactantExt/Simulations/simulation.jl +++ b/ext/OceananigansReactantExt/Simulations/simulation.jl @@ -31,4 +31,3 @@ time(sim::ReactantSimulation) = Reactant.to_number(time(sim.model)) add_callback!(::ReactantSimulation, args...) = error("Cannot add callbacks to a Simulation with ReactantState architecture!") - diff --git a/ext/OceananigansReactantExt/TimeSteppers.jl b/ext/OceananigansReactantExt/TimeSteppers.jl index 635d52fd985..b48b50e4447 100644 --- a/ext/OceananigansReactantExt/TimeSteppers.jl +++ b/ext/OceananigansReactantExt/TimeSteppers.jl @@ -134,4 +134,3 @@ function first_time_step!(model::ReactantModel{<:QuasiAdamsBashforth2TimeStepper end end # module - diff --git a/ext/OceananigansReactantExt/TurbulenceClosures.jl b/ext/OceananigansReactantExt/TurbulenceClosures.jl index 2670e8744ae..74e0f46e71b 100644 --- a/ext/OceananigansReactantExt/TurbulenceClosures.jl +++ b/ext/OceananigansReactantExt/TurbulenceClosures.jl @@ -8,4 +8,3 @@ import Oceananigans.TurbulenceClosures.TKEBasedVerticalDiffusivities: update_pre update_previous_compute_time!(diffusivities, model::ReactantModel) = model.clock.last_Δt end # module TurbulenceClosures - diff --git a/src/AbstractOperations/at.jl b/src/AbstractOperations/at.jl index f8fe2ecfc94..92c80ffbd93 100644 --- a/src/AbstractOperations/at.jl +++ b/src/AbstractOperations/at.jl @@ -27,7 +27,7 @@ end function instantiate_location_expression(exp::Expr) new_exp = deepcopy(exp) for (i, arg) in enumerate(new_exp.args) - new_exp.args[i] = instantiate_location_expression(arg) + new_exp.args[i] = instantiate_location_expression(arg) end return new_exp end diff --git a/src/AbstractOperations/binary_operations.jl b/src/AbstractOperations/binary_operations.jl index 9a4565ca480..df66e5c4e5b 100644 --- a/src/AbstractOperations/binary_operations.jl +++ b/src/AbstractOperations/binary_operations.jl @@ -126,7 +126,7 @@ function define_binary_operator(op) # instantiate location if types are passed $op(Lc::Tuple, a, b) = $op((Lc[1](), Lc[2](), Lc[3]()), a, b) - + $op(Lc::Tuple, f::Function, b::AbstractField) = $op((Lc[1](), Lc[2](), Lc[3]()), a, b) $op(Lc::Tuple, a::AbstractField, f::Function) = $op((Lc[1](), Lc[2](), Lc[3]()), a, b) @@ -186,8 +186,8 @@ BinaryOperation at (Center, Center, Center) ├── grid: 1×1×1 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 1×1×1 halo └── tree: plus_or_times at (Center, Center, Center) -    ├── 1×1×1 Field{Center, Center, Center} on RectilinearGrid on CPU -    └── 1×1×1 Field{Center, Center, Center} on RectilinearGrid on CPU + ├── 1×1×1 Field{Center, Center, Center} on RectilinearGrid on CPU + └── 1×1×1 Field{Center, Center, Center} on RectilinearGrid on CPU ``` """ macro binary(ops...) diff --git a/src/AbstractOperations/derivatives.jl b/src/AbstractOperations/derivatives.jl index 1b2b15a01e0..7b6eb43cf6a 100644 --- a/src/AbstractOperations/derivatives.jl +++ b/src/AbstractOperations/derivatives.jl @@ -85,8 +85,8 @@ by interpolation to `L`, where `L` is a 3-tuple of instantiated `Face`s and `Ce # Instantiate location if types are passed ∂x(L::Tuple, arg::AF) = ∂x((L[1](), L[2](), L[3]()), arg) -∂y(L::Tuple, arg::AF) = ∂y((L[1](), L[2](), L[3]()), arg) -∂z(L::Tuple, arg::AF) = ∂z((L[1](), L[2](), L[3]()), arg) +∂y(L::Tuple, arg::AF) = ∂y((L[1](), L[2](), L[3]()), arg) +∂z(L::Tuple, arg::AF) = ∂z((L[1](), L[2](), L[3]()), arg) # Defaults diff --git a/src/AbstractOperations/grid_metrics.jl b/src/AbstractOperations/grid_metrics.jl index 77ff27cf009..1e2e44a77d8 100644 --- a/src/AbstractOperations/grid_metrics.jl +++ b/src/AbstractOperations/grid_metrics.jl @@ -85,7 +85,7 @@ grid_metric_operation(Loc::Tuple, metric, grid) = grid_metric_operation((Loc[1]( const NodeMetric = Union{XNode, YNode, ZNode, ΛNode, ΦNode, RNode} -function grid_metric_operation(loc::Tuple{LX, LY, LZ}, metric::NodeMetric, grid) where {LX<:Location, LY<:Location, LZ<:Location} +function grid_metric_operation(loc::Tuple{LX, LY, LZ}, metric::NodeMetric, grid) where {LX<:Location, LY<:Location, LZ<:Location} ℓx, ℓy, ℓz = loc ξnode = metric_function(loc, metric) return KernelFunctionOperation{LX, LY, LZ}(ξnode, grid, ℓx, ℓy, ℓz) diff --git a/src/AbstractOperations/metric_field_reductions.jl b/src/AbstractOperations/metric_field_reductions.jl index 2af661c1ec8..2fe3be54b10 100644 --- a/src/AbstractOperations/metric_field_reductions.jl +++ b/src/AbstractOperations/metric_field_reductions.jl @@ -191,14 +191,14 @@ set!(c, z -> z) """ function CumulativeIntegral(field::AbstractField; dims, reverse=false, condition=nothing, mask=0) dims ∈ (1, 2, 3) || throw(ArgumentError("CumulativeIntegral only supports dims=1, 2, or 3.")) - + # Check that we're not accumulating over a Periodic dimension topo = topology(field.grid, dims) if topo === Periodic throw(ArgumentError("CumulativeIntegral does not support Periodic dimensions. " * "The cumulative integral of a periodic function is not periodic.")) end - + maybe_reverse_cumsum = reverse ? reverse_cumsum! : cumsum! dx = reduction_grid_metric(dims) operand = condition_operand(field * dx, condition, mask) diff --git a/src/AbstractOperations/multiary_operations.jl b/src/AbstractOperations/multiary_operations.jl index d949f81d1fa..788118636f8 100644 --- a/src/AbstractOperations/multiary_operations.jl +++ b/src/AbstractOperations/multiary_operations.jl @@ -22,33 +22,33 @@ end Π.op(Π.▶[1](i, j, k, Π.grid, Π.args[1])) @inline Base.getindex(Π::MultiaryOperation{LX, LY, LZ, 2}, i, j, k) where {LX, LY, LZ} = - Π.op(Π.▶[1](i, j, k, Π.grid, Π.args[1]), + Π.op(Π.▶[1](i, j, k, Π.grid, Π.args[1]), Π.▶[2](i, j, k, Π.grid, Π.args[2])) @inline Base.getindex(Π::MultiaryOperation{LX, LY, LZ, 3}, i, j, k) where {LX, LY, LZ} = - Π.op(Π.▶[1](i, j, k, Π.grid, Π.args[1]), - Π.▶[2](i, j, k, Π.grid, Π.args[2]), + Π.op(Π.▶[1](i, j, k, Π.grid, Π.args[1]), + Π.▶[2](i, j, k, Π.grid, Π.args[2]), Π.▶[3](i, j, k, Π.grid, Π.args[3])) @inline Base.getindex(Π::MultiaryOperation{LX, LY, LZ, 4}, i, j, k) where {LX, LY, LZ} = - Π.op(Π.▶[1](i, j, k, Π.grid, Π.args[1]), - Π.▶[2](i, j, k, Π.grid, Π.args[2]), - Π.▶[3](i, j, k, Π.grid, Π.args[3]), + Π.op(Π.▶[1](i, j, k, Π.grid, Π.args[1]), + Π.▶[2](i, j, k, Π.grid, Π.args[2]), + Π.▶[3](i, j, k, Π.grid, Π.args[3]), Π.▶[4](i, j, k, Π.grid, Π.args[4])) @inline Base.getindex(Π::MultiaryOperation{LX, LY, LZ, 5}, i, j, k) where {LX, LY, LZ} = - Π.op(Π.▶[1](i, j, k, Π.grid, Π.args[1]), - Π.▶[2](i, j, k, Π.grid, Π.args[2]), - Π.▶[3](i, j, k, Π.grid, Π.args[3]), - Π.▶[4](i, j, k, Π.grid, Π.args[4]), + Π.op(Π.▶[1](i, j, k, Π.grid, Π.args[1]), + Π.▶[2](i, j, k, Π.grid, Π.args[2]), + Π.▶[3](i, j, k, Π.grid, Π.args[3]), + Π.▶[4](i, j, k, Π.grid, Π.args[4]), Π.▶[5](i, j, k, Π.grid, Π.args[5])) @inline Base.getindex(Π::MultiaryOperation{LX, LY, LZ, 6}, i, j, k) where {LX, LY, LZ} = - Π.op(Π.▶[1](i, j, k, Π.grid, Π.args[1]), - Π.▶[2](i, j, k, Π.grid, Π.args[2]), - Π.▶[3](i, j, k, Π.grid, Π.args[3]), - Π.▶[4](i, j, k, Π.grid, Π.args[4]), - Π.▶[5](i, j, k, Π.grid, Π.args[5]), + Π.op(Π.▶[1](i, j, k, Π.grid, Π.args[1]), + Π.▶[2](i, j, k, Π.grid, Π.args[2]), + Π.▶[3](i, j, k, Π.grid, Π.args[3]), + Π.▶[4](i, j, k, Π.grid, Π.args[4]), + Π.▶[5](i, j, k, Π.grid, Π.args[5]), Π.▶[6](i, j, k, Π.grid, Π.args[6])) ##### @@ -85,10 +85,10 @@ function define_multiary_operator(op) end # Instantiate location if types are passed - $op(Lop::Tuple, - a::Union{Function, Number, Oceananigans.Fields.AbstractField}, - b::Union{Function, Number, Oceananigans.Fields.AbstractField}, - c::Union{Function, Number, Oceananigans.Fields.AbstractField}, + $op(Lop::Tuple, + a::Union{Function, Number, Oceananigans.Fields.AbstractField}, + b::Union{Function, Number, Oceananigans.Fields.AbstractField}, + c::Union{Function, Number, Oceananigans.Fields.AbstractField}, d::Union{Function, Number, Oceananigans.Fields.AbstractField}...) = $op((Lop[1](), Lop[2](), Lop[3]()), a, b, c, d...) $op(a::Oceananigans.Fields.AbstractField, @@ -125,17 +125,17 @@ BinaryOperation at (Center, Center, Center) ├── grid: 1×1×1 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 1×1×1 halo └── tree: * at (Center, Center, Center) -    ├── 0.3333333333333333 -    └── + at (Center, Center, Center) -       ├── / at (Center, Center, Center) -       │   ├── 1 -       │   └── 1×1×1 Field{Center, Center, Center} on RectilinearGrid on CPU -       ├── / at (Center, Center, Center) -       │   ├── 1 -       │   └── 1×1×1 Field{Center, Center, Center} on RectilinearGrid on CPU -       └── / at (Center, Center, Center) -          ├── 1 -          └── 1×1×1 Field{Center, Center, Center} on RectilinearGrid on CPU + ├── 0.3333333333333333 + └── + at (Center, Center, Center) + ├── / at (Center, Center, Center) + │ ├── 1 + │ └── 1×1×1 Field{Center, Center, Center} on RectilinearGrid on CPU + ├── / at (Center, Center, Center) + │ ├── 1 + │ └── 1×1×1 Field{Center, Center, Center} on RectilinearGrid on CPU + └── / at (Center, Center, Center) + ├── 1 + └── 1×1×1 Field{Center, Center, Center} on RectilinearGrid on CPU julia> @multiary harmonic_plus Set{Any} with 3 elements: @@ -148,9 +148,9 @@ MultiaryOperation at (Center, Center, Center) ├── grid: 1×1×1 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 1×1×1 halo └── tree: harmonic_plus at (Center, Center, Center) -    ├── 1×1×1 Field{Center, Center, Center} on RectilinearGrid on CPU -    ├── 1×1×1 Field{Center, Center, Center} on RectilinearGrid on CPU -    └── 1×1×1 Field{Center, Center, Center} on RectilinearGrid on CPU + ├── 1×1×1 Field{Center, Center, Center} on RectilinearGrid on CPU + ├── 1×1×1 Field{Center, Center, Center} on RectilinearGrid on CPU + └── 1×1×1 Field{Center, Center, Center} on RectilinearGrid on CPU ``` """ macro multiary(ops...) diff --git a/src/AbstractOperations/show_abstract_operations.jl b/src/AbstractOperations/show_abstract_operations.jl index 16f2f891024..a576ff92e18 100644 --- a/src/AbstractOperations/show_abstract_operations.jl +++ b/src/AbstractOperations/show_abstract_operations.jl @@ -32,7 +32,7 @@ tree_show(a::Union{Number, Function}, depth, nesting) = string(a) tree_show(a, depth, nesting) = summary(a) # fallback "Returns a string corresponding to padding characters for a tree visualization of an `AbstractOperation`." -get_tree_padding(depth, nesting) = "    "^(depth-nesting) * "│   "^nesting +get_tree_padding(depth, nesting) = " "^(depth-nesting) * "│ "^nesting "Return a string representaion of a `UnaryOperation` leaf within a tree visualization of an `AbstractOperation`." function tree_show(unary::UnaryOperation, depth, nesting) diff --git a/src/AbstractOperations/unary_operations.jl b/src/AbstractOperations/unary_operations.jl index c8279b24388..420c38c01cc 100644 --- a/src/AbstractOperations/unary_operations.jl +++ b/src/AbstractOperations/unary_operations.jl @@ -83,7 +83,7 @@ UnaryOperation at (Center, Center, Center) ├── grid: 1×1×1 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 1×1×1 halo └── tree: square_it at (Center, Center, Center) via identity -    └── 1×1×1 Field{Center, Center, Center} on RectilinearGrid on CPU + └── 1×1×1 Field{Center, Center, Center} on RectilinearGrid on CPU ``` """ macro unary(ops...) diff --git a/src/Advection/adapt_advection_order.jl b/src/Advection/adapt_advection_order.jl index d0b7f406330..1f2b98796c0 100644 --- a/src/Advection/adapt_advection_order.jl +++ b/src/Advection/adapt_advection_order.jl @@ -93,4 +93,4 @@ function adapt_advection_order(advection::WENO{B}, N::Int, grid::AbstractGrid) w else return WENO(order=2N-1) end -end \ No newline at end of file +end diff --git a/src/Advection/centered_reconstruction.jl b/src/Advection/centered_reconstruction.jl index df592fa7123..a2766742554 100644 --- a/src/Advection/centered_reconstruction.jl +++ b/src/Advection/centered_reconstruction.jl @@ -8,7 +8,7 @@ struct Centered{N, FT, CA} <: AbstractCenteredAdvectionScheme{N, FT} Centered{N, FT}(buffer_scheme::CA) where {N, FT, CA} = new{N, FT, CA}(buffer_scheme) end -function Centered(FT::DataType=Oceananigans.defaults.FloatType; +function Centered(FT::DataType=Oceananigans.defaults.FloatType; order = 2, buffer_scheme = DecreasingOrderAdvectionScheme()) diff --git a/src/Advection/flat_advective_fluxes.jl b/src/Advection/flat_advective_fluxes.jl index 9ada75ba96a..690ac88f0c7 100644 --- a/src/Advection/flat_advective_fluxes.jl +++ b/src/Advection/flat_advective_fluxes.jl @@ -49,4 +49,4 @@ for (dir, GridType) in zip((:xᶠᵃᵃ, :yᵃᶠᵃ, :zᵃᵃᶠ, :xᶜᵃᵃ, @inline $alt_biased_interp(i, j, k, grid::$GridType, ::HOADV, bias, ψ::Callable, ::AS, args...) = ψ(i, j, k, grid, args...) @inline $alt_biased_interp(i, j, k, grid::$GridType, ::LOADV, bias, ψ::Callable, ::AS, args...) = ψ(i, j, k, grid, args...) end -end \ No newline at end of file +end diff --git a/src/Advection/immersed_advective_fluxes.jl b/src/Advection/immersed_advective_fluxes.jl index 7c92c422997..35bb04e9c0c 100644 --- a/src/Advection/immersed_advective_fluxes.jl +++ b/src/Advection/immersed_advective_fluxes.jl @@ -188,7 +188,7 @@ for bias in (:symmetric, :biased) @inline $alt1_interp(i, j, k, ibg::ImmersedBoundaryGrid, scheme::LOADV, args...) = $interp(i, j, k, ibg, scheme, args...) # Conditional high-order interpolation in Bounded directions - @inline $alt1_interp(i, j, k, ibg::ImmersedBoundaryGrid, scheme::HOADV, args...) = + @inline $alt1_interp(i, j, k, ibg::ImmersedBoundaryGrid, scheme::HOADV, args...) = ifelse($near_boundary(i, j, k, ibg, scheme), $alt2_interp(i, j, k, ibg, scheme.buffer_scheme, args...), $interp(i, j, k, ibg, scheme, args...)) diff --git a/src/Advection/momentum_advection_operators.jl b/src/Advection/momentum_advection_operators.jl index fa334755b6f..fb2e5819680 100644 --- a/src/Advection/momentum_advection_operators.jl +++ b/src/Advection/momentum_advection_operators.jl @@ -97,4 +97,4 @@ end @inline _advective_momentum_flux_Wu(i, j, k, grid, ::Nothing, args...) = zero(grid) @inline _advective_momentum_flux_Wv(i, j, k, grid, ::Nothing, args...) = zero(grid) -@inline _advective_momentum_flux_Ww(i, j, k, grid, ::Nothing, args...) = zero(grid) \ No newline at end of file +@inline _advective_momentum_flux_Ww(i, j, k, grid, ::Nothing, args...) = zero(grid) diff --git a/src/Advection/multi_dimensional_reconstruction.jl b/src/Advection/multi_dimensional_reconstruction.jl index e57bb30202b..a610d80447b 100644 --- a/src/Advection/multi_dimensional_reconstruction.jl +++ b/src/Advection/multi_dimensional_reconstruction.jl @@ -131,4 +131,4 @@ end S₂ = (Q₀ , Q₊₁, Q₊₂) return fifth_order_weno_reconstruction(FT, S₀, S₁, S₂) -end \ No newline at end of file +end diff --git a/src/Advection/upwind_biased_reconstruction.jl b/src/Advection/upwind_biased_reconstruction.jl index 327934cfc0d..5c56e69991d 100644 --- a/src/Advection/upwind_biased_reconstruction.jl +++ b/src/Advection/upwind_biased_reconstruction.jl @@ -10,7 +10,7 @@ struct UpwindBiased{N, FT, CA, SI} <: AbstractUpwindBiasedAdvectionScheme{N, FT} new{N, FT, CA, SI}(buffer_scheme, advecting_velocity_scheme) end -function UpwindBiased(FT::DataType = Float64; +function UpwindBiased(FT::DataType = Float64; order = 3, buffer_scheme = DecreasingOrderAdvectionScheme(), minimum_buffer_upwind_order = 1) @@ -48,7 +48,7 @@ Base.summary(a::UpwindBiased{N}) where N = string("UpwindBiased(order=", 2N-1, " function print_buffer_scheme_tree(io::IO, scheme, prefix::String, is_last::Bool) connector = is_last ? "└── " : "├── " print(io, prefix, connector, "buffer_scheme: ", summary(scheme)) - + # Check if this scheme has a nested buffer_scheme to display if hasproperty(scheme, :buffer_scheme) && !isnothing(scheme.buffer_scheme) println(io) diff --git a/src/Advection/weno_reconstruction.jl b/src/Advection/weno_reconstruction.jl index 4f98b1e7f79..65a90346923 100644 --- a/src/Advection/weno_reconstruction.jl +++ b/src/Advection/weno_reconstruction.jl @@ -154,4 +154,3 @@ on_architecture(to, scheme::WENO{N, FT, FT2}) where {N, FT, FT2} = WENO{N, FT, FT2}(on_architecture(to, scheme.bounds), on_architecture(to, scheme.buffer_scheme), on_architecture(to, scheme.advecting_velocity_scheme)) - diff --git a/src/BoundaryConditions/discrete_boundary_function.jl b/src/BoundaryConditions/discrete_boundary_function.jl index 57bda13d662..fbf3f8bab19 100644 --- a/src/BoundaryConditions/discrete_boundary_function.jl +++ b/src/BoundaryConditions/discrete_boundary_function.jl @@ -57,4 +57,3 @@ Adapt.adapt_structure(to, bf::DiscreteBoundaryFunction) = DiscreteBoundaryFuncti on_architecture(to, bf::DiscreteBoundaryFunction) = DiscreteBoundaryFunction(on_architecture(to, bf.func), on_architecture(to, bf.parameters)) - diff --git a/src/BoundaryConditions/fill_halo_regions_periodic.jl b/src/BoundaryConditions/fill_halo_regions_periodic.jl index 4c51a493576..9ca3dc3714a 100644 --- a/src/BoundaryConditions/fill_halo_regions_periodic.jl +++ b/src/BoundaryConditions/fill_halo_regions_periodic.jl @@ -30,4 +30,4 @@ end parent(c)[i, j, k] = parent(c)[i, j, N+k] # top parent(c)[i, j, N+H+k] = parent(c)[i, j, H+k] # bottom end -end \ No newline at end of file +end diff --git a/src/BoundaryConditions/polar_boundary_condition.jl b/src/BoundaryConditions/polar_boundary_condition.jl index a74bee285fc..a80825f3f39 100644 --- a/src/BoundaryConditions/polar_boundary_condition.jl +++ b/src/BoundaryConditions/polar_boundary_condition.jl @@ -79,4 +79,3 @@ function fill_halo_event!(c, kernel!, bcs::SouthAndNorthPolarBC, loc, grid, args update_pole_value!(bcs[2].condition, c, grid, loc) return kernel!(c, bcs..., loc, grid, Tuple(args)) end - diff --git a/src/BuoyancyFormulations/buoyancy_force.jl b/src/BuoyancyFormulations/buoyancy_force.jl index c5542a43e85..cce42d2013f 100644 --- a/src/BuoyancyFormulations/buoyancy_force.jl +++ b/src/BuoyancyFormulations/buoyancy_force.jl @@ -66,12 +66,12 @@ function BuoyancyForce(grid, formulation::AbstractBuoyancyFormulation; gravity_u end # Fallback for when no grid is available, we overwrite `materialize_gradients` to false -BuoyancyForce(formulation::AbstractBuoyancyFormulation; materialize_gradients=false, kwargs...) = +BuoyancyForce(formulation::AbstractBuoyancyFormulation; materialize_gradients=false, kwargs...) = BuoyancyForce(nothing, formulation; materialize_gradients=false, kwargs...) -Adapt.adapt_structure(to, bf::BuoyancyForce) = - BuoyancyForce(Adapt.adapt(to, bf.formulation), - Adapt.adapt(to, bf.gravity_unit_vector), +Adapt.adapt_structure(to, bf::BuoyancyForce) = + BuoyancyForce(Adapt.adapt(to, bf.formulation), + Adapt.adapt(to, bf.gravity_unit_vector), Adapt.adapt(to, bf.gradients)) @inline ĝ_x(bf) = @inbounds - bf.gravity_unit_vector[1] @@ -101,7 +101,7 @@ materialize_buoyancy(formulation::AbstractBuoyancyFormulation, grid; kw...) = Bu # Fallback compute_buoyancy_gradients!(::BuoyancyForce{<:Any, <:Any, <:Nothing}, grid, tracers; kw...) = nothing -compute_buoyancy_gradients!(::Nothing, grid, tracers; kw...) = nothing +compute_buoyancy_gradients!(::Nothing, grid, tracers; kw...) = nothing Base.summary(bf::BuoyancyForce) = string(summary(bf.formulation), " with ĝ = ", @@ -127,7 +127,7 @@ end @inline ∂y_b(i, j, k, grid, b::BuoyancyForce, C) = @inbounds b.gradients.∂y_b[i, j, k] @inline ∂z_b(i, j, k, grid, b::BuoyancyForce, C) = @inbounds b.gradients.∂z_b[i, j, k] -function compute_buoyancy_gradients!(buoyancy, grid, tracers; parameters=:xyz) +function compute_buoyancy_gradients!(buoyancy, grid, tracers; parameters=:xyz) gradients = buoyancy.gradients formulation = buoyancy.formulation launch!(architecture(grid), grid, parameters, _compute_buoyancy_gradients!, gradients, grid, formulation, tracers) diff --git a/src/BuoyancyFormulations/g_dot_b.jl b/src/BuoyancyFormulations/g_dot_b.jl index c3a4eaa6977..ded5c1b910f 100644 --- a/src/BuoyancyFormulations/g_dot_b.jl +++ b/src/BuoyancyFormulations/g_dot_b.jl @@ -8,4 +8,3 @@ const NZBF = BuoyancyForce{<:Any, NegativeZDirection} @inline x_dot_g_bᶠᶜᶜ(i, j, k, grid, bf::NZBF, C) = 0 @inline y_dot_g_bᶜᶠᶜ(i, j, k, grid, bf::NZBF, C) = 0 @inline z_dot_g_bᶜᶜᶠ(i, j, k, grid, bf::NZBF, C) = ℑzᵃᵃᶠ(i, j, k, grid, buoyancy_perturbationᶜᶜᶜ, bf.formulation, C) - diff --git a/src/BuoyancyFormulations/linear_equation_of_state.jl b/src/BuoyancyFormulations/linear_equation_of_state.jl index 23d17f4f846..af60bcabc84 100644 --- a/src/BuoyancyFormulations/linear_equation_of_state.jl +++ b/src/BuoyancyFormulations/linear_equation_of_state.jl @@ -78,4 +78,3 @@ const LinearSalinitySeawaterBuoyancy = SeawaterBuoyancy{FT, <:LinearEquationOfSt @inline buoyancy_perturbationᶜᶜᶜ(i, j, k, grid, b::LinearSalinitySeawaterBuoyancy, C) = @inbounds - b.gravitational_acceleration * b.equation_of_state.haline_contraction * C.S[i, j, k] - diff --git a/src/BuoyancyFormulations/seawater_buoyancy.jl b/src/BuoyancyFormulations/seawater_buoyancy.jl index 44359fa2083..0f1f3a53e05 100644 --- a/src/BuoyancyFormulations/seawater_buoyancy.jl +++ b/src/BuoyancyFormulations/seawater_buoyancy.jl @@ -245,4 +245,3 @@ end @inline top_buoyancy_flux(i, j, grid, b::SeawaterBuoyancy, args...) = top_bottom_buoyancy_flux(i, j, grid.Nz+1, grid, b, args...) @inline bottom_buoyancy_flux(i, j, grid, b::SeawaterBuoyancy, args...) = top_bottom_buoyancy_flux(i, j, 1, grid, b, args...) - diff --git a/src/Coriolis/f_plane.jl b/src/Coriolis/f_plane.jl index 6e588d58052..166c31a4c42 100644 --- a/src/Coriolis/f_plane.jl +++ b/src/Coriolis/f_plane.jl @@ -57,4 +57,3 @@ function Base.summary(fplane::FPlane{FT}) where FT end Base.show(io::IO, fplane::FPlane) = print(io, summary(fplane)) - diff --git a/src/DistributedComputations/communication_buffers.jl b/src/DistributedComputations/communication_buffers.jl index 1c9862406d5..7b4b354ed59 100644 --- a/src/DistributedComputations/communication_buffers.jl +++ b/src/DistributedComputations/communication_buffers.jl @@ -26,23 +26,23 @@ struct CommunicationBuffers{W, E, S, N, SW, SE, NW, NE} end Adapt.adapt_structure(to, buff::CommunicationBuffers) = - CommunicationBuffers(Adapt.adapt(to, buff.west), - Adapt.adapt(to, buff.east), - Adapt.adapt(to, buff.north), - Adapt.adapt(to, buff.south), - Adapt.adapt(to, buff.southwest), - Adapt.adapt(to, buff.southeast), - Adapt.adapt(to, buff.northwest), + CommunicationBuffers(Adapt.adapt(to, buff.west), + Adapt.adapt(to, buff.east), + Adapt.adapt(to, buff.north), + Adapt.adapt(to, buff.south), + Adapt.adapt(to, buff.southwest), + Adapt.adapt(to, buff.southeast), + Adapt.adapt(to, buff.northwest), Adapt.adapt(to, buff.northeast)) on_architecture(arch, buff::CommunicationBuffers) = - CommunicationBuffers(on_architecture(arch, buff.west), - on_architecture(arch, buff.east), - on_architecture(arch, buff.north), - on_architecture(arch, buff.south), - on_architecture(arch, buff.southwest), - on_architecture(arch, buff.southeast), - on_architecture(arch, buff.northwest), + CommunicationBuffers(on_architecture(arch, buff.west), + on_architecture(arch, buff.east), + on_architecture(arch, buff.north), + on_architecture(arch, buff.south), + on_architecture(arch, buff.southwest), + on_architecture(arch, buff.southeast), + on_architecture(arch, buff.northwest), on_architecture(arch, buff.northeast)) communication_buffers(grid::DistributedGrid, data, boundary_conditions) = CommunicationBuffers(grid, data, boundary_conditions) @@ -52,9 +52,9 @@ communication_buffers(grid::DistributedGrid, data, boundary_conditions) = Commun Construct communication buffers for distributed halo exchange. -`CommunicationBuffers` stores send/receive buffers for each spatial direction and corner -in a distributed grid. During halo exchange, data is copied from the interior domain into -send buffers, communicated via MPI to neighboring processes, and then unpacked from receive +`CommunicationBuffers` stores send/receive buffers for each spatial direction and corner +in a distributed grid. During halo exchange, data is copied from the interior domain into +send buffers, communicated via MPI to neighboring processes, and then unpacked from receive buffers into halo regions. # Buffer Types @@ -96,12 +96,12 @@ CommunicationBuffers(grid, data, ::Nothing) = nothing Communication buffer for one-dimensional domain decomposition or edge boundaries. -In a one-dimensional parallelization (e.g., only in x or only in y), `OneDBuffer` -contains the full extent of the halo region including the corners. This allows corner +In a one-dimensional parallelization (e.g., only in x or only in y), `OneDBuffer` +contains the full extent of the halo region including the corners. This allows corner data to be communicated in a single pass along with the edge data. -`OneDBuffer` is also used in two-dimensional parallelizations for processes at domain -edges (e.g., boundaries with `RightConnected` or `LeftConnected` topologies), where +`OneDBuffer` is also used in two-dimensional parallelizations for processes at domain +edges (e.g., boundaries with `RightConnected` or `LeftConnected` topologies), where corner communication is not needed in that direction. # Size @@ -118,9 +118,9 @@ end Communication buffer for two-dimensional domain decomposition without corner regions. -In a two-dimensional parallelization where corners are communicated separately via -`CornerBuffer`, `TwoDBuffer` contains only the edge halo region excluding the corners. -This enables efficient communication patterns where edge and corner data are sent in +In a two-dimensional parallelization where corners are communicated separately via +`CornerBuffer`, `TwoDBuffer` contains only the edge halo region excluding the corners. +This enables efficient communication patterns where edge and corner data are sent in separate MPI passes. # Size @@ -137,17 +137,17 @@ end Communication buffer for corner regions in two-dimensional domain decomposition. -In a two-dimensional parallelization, `CornerBuffer` handles the communication of -diagonal corner regions (southwest, southeast, northwest, northeast) that are not -covered by the edge communication buffers (`TwoDBuffer`). Corner communication ensures -that all halo regions are properly synchronized between neighboring processes in both +In a two-dimensional parallelization, `CornerBuffer` handles the communication of +diagonal corner regions (southwest, southeast, northwest, northeast) that are not +covered by the edge communication buffers (`TwoDBuffer`). Corner communication ensures +that all halo regions are properly synchronized between neighboring processes in both x and y directions. # Size `(Hx, Hy, Tz)` where `Hx` is the halo size in x, `Hy` is the halo size in y, and `Tz` is the size in z # Note -Corner buffers are only created for `Distributed` architectures with two-dimensional +Corner buffers are only created for `Distributed` architectures with two-dimensional parallelization and are `nothing` otherwise. """ struct CornerBuffer{B} @@ -199,7 +199,7 @@ function y_communication_buffer(arch::Distributed, grid::AbstractGrid{<:Any, TX, end # Never pass corners in a MCBC. -function x_communication_buffer(arch, grid, data, H, ::MCBC) +function x_communication_buffer(arch, grid, data, H, ::MCBC) _, Ty, Tz = size(parent(data)) FT = eltype(data) send = on_architecture(arch, zeros(FT, H, Ty, Tz)) @@ -207,7 +207,7 @@ function x_communication_buffer(arch, grid, data, H, ::MCBC) return OneDBuffer(send, recv) end -function y_communication_buffer(arch, grid, data, H, ::MCBC) +function y_communication_buffer(arch, grid, data, H, ::MCBC) Tx, _, Tz = size(parent(data)) FT = eltype(data) send = on_architecture(arch, zeros(FT, Tx, H, Tz)) @@ -226,7 +226,7 @@ corner_communication_buffer(arch, grid, data, Hx, Hy, xedge, yedge) = nothing corner_communication_buffer(::Distributed, grid, data, Hx, Hy, ::Nothing, ::Nothing) = nothing corner_communication_buffer(arch::Distributed, grid, data, Hx, Hy, ::Nothing, yedge) = nothing corner_communication_buffer(arch::Distributed, grid, data, Hx, Hy, xedge, ::Nothing) = nothing - + # CornerBuffer are used only in the two-dimensional partitioning case, in all other cases they are equal to `nothing` function corner_communication_buffer(arch::Distributed, grid, data, Hx, Hy, xedge, yedge) Tz = size(parent(data), 3) @@ -277,13 +277,13 @@ end ##### Single sided fill_send_buffers! ##### -fill_send_buffers!(c::OffsetArray, buff::CommunicationBuffers, grid, ::West) = +fill_send_buffers!(c::OffsetArray, buff::CommunicationBuffers, grid, ::West) = _fill_west_send_buffer!(parent(c), buff.west, halo_size(grid)[[1, 2]]..., size(grid)[[1, 2]]...) -fill_send_buffers!(c::OffsetArray, buff::CommunicationBuffers, grid, ::East) = +fill_send_buffers!(c::OffsetArray, buff::CommunicationBuffers, grid, ::East) = _fill_east_send_buffer!(parent(c), buff.east, halo_size(grid)[[1, 2]]..., size(grid)[[1, 2]]...) -fill_send_buffers!(c::OffsetArray, buff::CommunicationBuffers, grid, ::South) = +fill_send_buffers!(c::OffsetArray, buff::CommunicationBuffers, grid, ::South) = _fill_south_send_buffer!(parent(c), buff.south, halo_size(grid)[[1, 2]]..., size(grid)[[1, 2]]...) -fill_send_buffers!(c::OffsetArray, buff::CommunicationBuffers, grid, ::North) = +fill_send_buffers!(c::OffsetArray, buff::CommunicationBuffers, grid, ::North) = _fill_north_send_buffer!(parent(c), buff.north, halo_size(grid)[[1, 2]]..., size(grid)[[1, 2]]...) fill_send_buffers!(c::OffsetArray, buff::CommunicationBuffers, grid, ::Bottom) = nothing fill_send_buffers!(c::OffsetArray, buff::CommunicationBuffers, grid, ::Top) = nothing @@ -348,13 +348,13 @@ end ##### Single sided recv_from_buffers! ##### -recv_from_buffers!(c::OffsetArray, buff::CommunicationBuffers, grid, ::West) = +recv_from_buffers!(c::OffsetArray, buff::CommunicationBuffers, grid, ::West) = _recv_from_west_buffer!(parent(c), buff.west, halo_size(grid)[[1, 2]]..., size(grid)[[1, 2]]...) -recv_from_buffers!(c::OffsetArray, buff::CommunicationBuffers, grid, ::East) = +recv_from_buffers!(c::OffsetArray, buff::CommunicationBuffers, grid, ::East) = _recv_from_east_buffer!(parent(c), buff.east, halo_size(grid)[[1, 2]]..., size(grid)[[1, 2]]...) -recv_from_buffers!(c::OffsetArray, buff::CommunicationBuffers, grid, ::South) = +recv_from_buffers!(c::OffsetArray, buff::CommunicationBuffers, grid, ::South) = _recv_from_south_buffer!(parent(c), buff.south, halo_size(grid)[[1, 2]]..., size(grid)[[1, 2]]...) -recv_from_buffers!(c::OffsetArray, buff::CommunicationBuffers, grid, ::North) = +recv_from_buffers!(c::OffsetArray, buff::CommunicationBuffers, grid, ::North) = _recv_from_north_buffer!(parent(c), buff.north, halo_size(grid)[[1, 2]]..., size(grid)[[1, 2]]...) recv_from_buffers!(c::OffsetArray, buff::CommunicationBuffers, grid, ::Bottom) = nothing recv_from_buffers!(c::OffsetArray, buff::CommunicationBuffers, grid, ::Top) = nothing diff --git a/src/DistributedComputations/distributed_fields.jl b/src/DistributedComputations/distributed_fields.jl index 1ee05defd83..e7c99f887fe 100644 --- a/src/DistributedComputations/distributed_fields.jl +++ b/src/DistributedComputations/distributed_fields.jl @@ -72,7 +72,7 @@ synchronize_communication!(var) = throw(ArgumentError("`synchronize_communicatio # Methods for types that do not require synchronization synchronize_communication!(::AbstractField) = nothing -synchronize_communication!(::AbstractArray) = nothing +synchronize_communication!(::AbstractArray) = nothing synchronize_communication!(::Number) = nothing synchronize_communication!(::Nothing) = nothing @@ -228,14 +228,14 @@ end # Distributed dot product @inline function dot(u::DistributedField, v::DistributedField; condition=nothing) - cu = condition_operand(u, condition, 0) - cv = condition_operand(v, condition, 0) - - B = cu * cv # Binary operation - r = zeros(u.grid, 1) - - Base.mapreducedim!(identity, +, r, B) - dot_local = @allowscalar r[1] + cu = condition_operand(u, condition, 0) + cv = condition_operand(v, condition, 0) + + B = cu * cv # Binary operation + r = zeros(u.grid, 1) + + Base.mapreducedim!(identity, +, r, B) + dot_local = @allowscalar r[1] arch = architecture(u) return all_reduce(+, dot_local, arch) end @@ -269,4 +269,4 @@ end @inline mean(f::Function, c::DistributedAbstractField, dims; condition=nothing, mask=0) = _mean(f, c, dims; condition, mask) -@inline mean(c::DistributedAbstractField; condition=nothing, dims=:) = _mean(identity, c, dims; condition) \ No newline at end of file +@inline mean(c::DistributedAbstractField; condition=nothing, dims=:) = _mean(identity, c, dims; condition) diff --git a/src/DistributedComputations/distributed_kernel_launching.jl b/src/DistributedComputations/distributed_kernel_launching.jl index 0f17f7ebae5..9ff66fbcf75 100644 --- a/src/DistributedComputations/distributed_kernel_launching.jl +++ b/src/DistributedComputations/distributed_kernel_launching.jl @@ -4,4 +4,3 @@ function _launch!(arch::Distributed, args...; kwargs...) child_arch = child_architecture(arch) return _launch!(child_arch, args...; kwargs...) end - diff --git a/src/DistributedComputations/distributed_on_architecture.jl b/src/DistributedComputations/distributed_on_architecture.jl index 00920a257eb..82c0491262f 100644 --- a/src/DistributedComputations/distributed_on_architecture.jl +++ b/src/DistributedComputations/distributed_on_architecture.jl @@ -74,4 +74,4 @@ function on_architecture(new_arch::Distributed, old_grid::OrthogonalSphericalShe horizontal_area_data..., old_grid.radius, old_grid.conformal_mapping) -end \ No newline at end of file +end diff --git a/src/DistributedComputations/plan_distributed_transforms.jl b/src/DistributedComputations/plan_distributed_transforms.jl index f7f4b96a3a0..8fb7e8d18db 100644 --- a/src/DistributedComputations/plan_distributed_transforms.jl +++ b/src/DistributedComputations/plan_distributed_transforms.jl @@ -39,4 +39,4 @@ function plan_distributed_transforms(global_grid, storage::TransposableField, pl ) return (; forward = forward_operations, backward = backward_operations) -end \ No newline at end of file +end diff --git a/src/Fields/Fields.jl b/src/Fields/Fields.jl index 3298274e209..585631b5f5d 100644 --- a/src/Fields/Fields.jl +++ b/src/Fields/Fields.jl @@ -57,7 +57,7 @@ Build a field from array `a` at `loc` and on `grid`. end # Build a field off of the current data -@inline function field(loc, a::OffsetArray, grid) +@inline function field(loc, a::OffsetArray, grid) loc = instantiate(loc) a = on_architecture(architecture(grid), a) return Field(loc, grid; data=a) diff --git a/src/Fields/constant_field.jl b/src/Fields/constant_field.jl index 3e05ea4135b..726e0503257 100644 --- a/src/Fields/constant_field.jl +++ b/src/Fields/constant_field.jl @@ -21,4 +21,3 @@ const CF = Union{ConstantField, ZeroField, OneField} BoundaryConditions.fill_halo_regions!(::ZeroField, args...; kw...) = nothing BoundaryConditions.fill_halo_regions!(::ConstantField, args...; kw...) = nothing - diff --git a/src/Fields/field_tuples.jl b/src/Fields/field_tuples.jl index f8884d356df..2974991f220 100644 --- a/src/Fields/field_tuples.jl +++ b/src/Fields/field_tuples.jl @@ -135,7 +135,7 @@ function VelocityFields(grid::AbstractGrid, user_bcs = NamedTuple()) u = XFaceField(grid, boundary_conditions=bcs.u) v = YFaceField(grid, boundary_conditions=bcs.v) w = ZFaceField(grid, boundary_conditions=bcs.w) - + return (u=u, v=v, w=w) end diff --git a/src/Fields/function_field.jl b/src/Fields/function_field.jl index 0d63fdf2aae..f06d76ef9de 100644 --- a/src/Fields/function_field.jl +++ b/src/Fields/function_field.jl @@ -78,4 +78,3 @@ Base.show(io::IO, field::FunctionField) = "├── grid: $(summary(field.grid))\n", "├── clock: $(summary(field.clock))\n", "└── parameters: $(field.parameters)") - diff --git a/src/Fields/interpolate.jl b/src/Fields/interpolate.jl index edf48b07b77..a5d1c541fe7 100644 --- a/src/Fields/interpolate.jl +++ b/src/Fields/interpolate.jl @@ -408,4 +408,3 @@ function interpolate!(to_field::Field, from_field::AbstractField) return to_field end - diff --git a/src/Fields/scans.jl b/src/Fields/scans.jl index 5a8bcf5b6eb..27fa2544e99 100644 --- a/src/Fields/scans.jl +++ b/src/Fields/scans.jl @@ -203,7 +203,7 @@ Accumulation(accumulate!, operand; dims) = Scan(Accumulating(), accumulate!, ope flip(::Type{Face}) = Center flip(::Type{Center}) = Face - + function Oceananigans.location(a::Accumulation) op_loc = location(a.operand) loc = Tuple(d ∈ a.dims ? flip(op_loc[d]) : op_loc[d] for d=1:3) diff --git a/src/Fields/set!.jl b/src/Fields/set!.jl index 81e160de84c..9d541da43a4 100644 --- a/src/Fields/set!.jl +++ b/src/Fields/set!.jl @@ -92,7 +92,7 @@ function set_to_function!(u, f, clock=nothing) # Transfer data to GPU if u is on the GPU if child_arch isa GPU || child_arch isa ReactantState - set!(u, cpu_u) + set!(u, cpu_u) end return u end @@ -151,4 +151,3 @@ end Base.copyto!(f::Field, src::Base.Broadcast.Broadcasted) = copyto!(interior(f), src) Base.copyto!(f::Field, src::AbstractArray) = copyto!(interior(f), src) Base.copyto!(f::Field, src::Field) = copyto!(parent(f), parent(src)) - diff --git a/src/Forcings/forcing.jl b/src/Forcings/forcing.jl index 57904871df8..8f5c66da58b 100644 --- a/src/Forcings/forcing.jl +++ b/src/Forcings/forcing.jl @@ -184,4 +184,3 @@ Return a `Forcing` by a `FieldTimeSeries`, which can be added to the tendency of Forcing is computed by calling `fts[i, j, k, Time(clock.time)]`, so the `FieldTimeSeries` must have the spatial dimensions of the `grid`. """ Forcing(fts::FlavorOfFTS) = Forcing(field_time_series_forcing_func; discrete_form=true, parameters=fts) - diff --git a/src/Forcings/model_forcing.jl b/src/Forcings/model_forcing.jl index aaadee52252..9896e37c7a9 100644 --- a/src/Forcings/model_forcing.jl +++ b/src/Forcings/model_forcing.jl @@ -47,7 +47,7 @@ function model_forcing(user_forcings, model_fields, prognostic_fields=model_fiel if name ∉ keys(prognostic_fields) msg = string("Invalid forcing: forcing contains an entry for $name, but $name is not a prognostic field!", '\n', "The prognostic fields are ", keys(prognostic_fields)) - throw(ArgumentError(msg)) + throw(ArgumentError(msg)) end end diff --git a/src/Forcings/multiple_forcings.jl b/src/Forcings/multiple_forcings.jl index 94dffcdb779..c3e2ed52ee6 100644 --- a/src/Forcings/multiple_forcings.jl +++ b/src/Forcings/multiple_forcings.jl @@ -72,4 +72,3 @@ function Base.show(io::IO, mf::MultipleForcings) return nothing end - diff --git a/src/Grids/latitude_longitude_grid.jl b/src/Grids/latitude_longitude_grid.jl index 2d3fbb5145f..c02708d4136 100644 --- a/src/Grids/latitude_longitude_grid.jl +++ b/src/Grids/latitude_longitude_grid.jl @@ -416,9 +416,9 @@ function Adapt.adapt_structure(to, grid::LatitudeLongitudeGrid) return LatitudeLongitudeGrid{TX, TY, TZ}(nothing, grid.Nx, grid.Ny, grid.Nz, grid.Hx, grid.Hy, grid.Hz, - Adapt.adapt(to, grid.Lx), - Adapt.adapt(to, grid.Ly), - Adapt.adapt(to, grid.Lz), + Adapt.adapt(to, grid.Lx), + Adapt.adapt(to, grid.Ly), + Adapt.adapt(to, grid.Lz), Adapt.adapt(to, grid.Δλᶠᵃᵃ), Adapt.adapt(to, grid.Δλᶜᵃᵃ), Adapt.adapt(to, grid.λᶠᵃᵃ), @@ -438,7 +438,7 @@ function Adapt.adapt_structure(to, grid::LatitudeLongitudeGrid) Adapt.adapt(to, grid.Azᶠᶜᵃ), Adapt.adapt(to, grid.Azᶜᶠᵃ), Adapt.adapt(to, grid.Azᶠᶠᵃ), - Adapt.adapt(to, grid.radius)) + Adapt.adapt(to, grid.radius)) end ##### diff --git a/src/Grids/new_data.jl b/src/Grids/new_data.jl index 70e29e474d6..3d7848490d6 100644 --- a/src/Grids/new_data.jl +++ b/src/Grids/new_data.jl @@ -74,4 +74,3 @@ new_data(FT::DataType, grid::AbstractGrid, loc, indices=default_indices(length(l new_data(FT, architecture(grid), loc, topology(grid), size(grid), halo_size(grid), indices) new_data(grid::AbstractGrid, loc, indices=default_indices) = new_data(eltype(grid), grid, loc, indices) - diff --git a/src/Grids/zeros_and_ones.jl b/src/Grids/zeros_and_ones.jl index 34529de2e68..a807edeeae7 100644 --- a/src/Grids/zeros_and_ones.jl +++ b/src/Grids/zeros_and_ones.jl @@ -10,4 +10,3 @@ Base.zeros(grid::AbstractGrid, N...) = zeros(architecture(grid), eltype(grid), N @inline Base.zero(grid::AbstractGrid) = zero(eltype(grid)) @inline Base.one(grid::AbstractGrid) = one(eltype(grid)) - diff --git a/src/ImmersedBoundaries/mutable_immersed_grid.jl b/src/ImmersedBoundaries/mutable_immersed_grid.jl index 7496c3ca34d..9fcd9310738 100644 --- a/src/ImmersedBoundaries/mutable_immersed_grid.jl +++ b/src/ImmersedBoundaries/mutable_immersed_grid.jl @@ -66,4 +66,4 @@ for LX in (:ᶠ, :ᶜ), LY in (:ᶠ, :ᶜ), LZ in (:ᶠ, :ᶜ) @inline $zspacing(i, j, k, grid::IMLLG) = $rspacing(i, j, k, grid) * σⁿ(i, j, k, grid, $ℓx(), $ℓy(), $ℓz()) @inline $zspacing(i, j, k, grid::IMOSG) = $rspacing(i, j, k, grid) * σⁿ(i, j, k, grid, $ℓx(), $ℓy(), $ℓz()) end -end \ No newline at end of file +end diff --git a/src/Models/HydrostaticFreeSurfaceModels/SplitExplicitFreeSurfaces/distributed_split_explicit_free_surface.jl b/src/Models/HydrostaticFreeSurfaceModels/SplitExplicitFreeSurfaces/distributed_split_explicit_free_surface.jl index cf8b067253a..b0d84ecfcf7 100644 --- a/src/Models/HydrostaticFreeSurfaceModels/SplitExplicitFreeSurfaces/distributed_split_explicit_free_surface.jl +++ b/src/Models/HydrostaticFreeSurfaceModels/SplitExplicitFreeSurfaces/distributed_split_explicit_free_surface.jl @@ -33,6 +33,6 @@ function synchronize_communication!(free_surface::SplitExplicitFreeSurface) for field in (U, V, Ũ, Ṽ, η) synchronize_communication!(field) end - + return nothing end diff --git a/src/Models/HydrostaticFreeSurfaceModels/SplitExplicitFreeSurfaces/step_split_explicit_free_surface.jl b/src/Models/HydrostaticFreeSurfaceModels/SplitExplicitFreeSurfaces/step_split_explicit_free_surface.jl index 55610f3ba69..a921abbdc29 100644 --- a/src/Models/HydrostaticFreeSurfaceModels/SplitExplicitFreeSurfaces/step_split_explicit_free_surface.jl +++ b/src/Models/HydrostaticFreeSurfaceModels/SplitExplicitFreeSurfaces/step_split_explicit_free_surface.jl @@ -13,12 +13,12 @@ Hᶠᶜ = column_depthᶠᶜᵃ(i, j, k_top, grid, η) Hᶜᶠ = column_depthᶜᶠᵃ(i, j, k_top, grid, η) - + # ∂τ(U) = - ∇η + G @inbounds begin U[i, j, 1] += Δτ * (- g * Hᶠᶜ * ∂xTᶠᶜᶠ(i, j, k_top, grid, η★, timestepper, η) + Gᵁ[i, j, 1]) V[i, j, 1] += Δτ * (- g * Hᶜᶠ * ∂yTᶜᶠᶠ(i, j, k_top, grid, η★, timestepper, η) + Gⱽ[i, j, 1]) - + # averaging the transport Ũ[i, j, 1] += transport_weight * U[i, j, 1] Ṽ[i, j, 1] += transport_weight * V[i, j, 1] @@ -32,7 +32,7 @@ end cache_previous_free_surface!(timestepper, i, j, k_top, η) δh_U = (δxTᶜᵃᵃ(i, j, grid.Nz, grid, Δy_qᶠᶜᶠ, U★, timestepper, U) + - δyTᵃᶜᵃ(i, j, grid.Nz, grid, Δx_qᶜᶠᶠ, U★, timestepper, V)) * Az⁻¹ᶜᶜᶠ(i, j, k_top, grid) + δyTᵃᶜᵃ(i, j, grid.Nz, grid, Δx_qᶜᶠᶠ, U★, timestepper, V)) * Az⁻¹ᶜᶜᶠ(i, j, k_top, grid) @inbounds begin η[i, j, k_top] += Δτ * (F(i, j, k_top, grid, clock, (; η, U, V)) - δh_U) @@ -91,7 +91,7 @@ function iterate_split_explicit!(free_surface, grid, GUⁿ, GVⁿ, Δτᴮ, F, c for substep in 1:Nsubsteps @inbounds averaging_weight = weights[substep] @inbounds transport_weight = transport_weights[substep] - + barotropic_velocity_kernel!(transport_weight, converted_U_args...) free_surface_kernel!(averaging_weight, converted_η_args...) end @@ -149,11 +149,11 @@ function step_free_surface!(free_surface::SplitExplicitFreeSurface, model, baroc # Wait for setup step to finish wait_free_surface_communication!(free_surface, model, architecture(free_surface_grid)) - + @apply_regionally begin - # Reset the filtered fields and the barotropic timestepper to zero. + # Reset the filtered fields and the barotropic timestepper to zero. initialize_free_surface_state!(free_surface, baroclinic_timestepper, barotropic_timestepper) - + # Solve for the free surface at tⁿ⁺¹ iterate_split_explicit!(free_surface, free_surface_grid, GUⁿ, GVⁿ, Δτᴮ, F, model.clock, weights, transport_weights, Val(Nsubsteps)) diff --git a/src/Models/HydrostaticFreeSurfaceModels/cache_hydrostatic_free_surface_tendencies.jl b/src/Models/HydrostaticFreeSurfaceModels/cache_hydrostatic_free_surface_tendencies.jl index 8ada686188a..5bb439ecc35 100644 --- a/src/Models/HydrostaticFreeSurfaceModels/cache_hydrostatic_free_surface_tendencies.jl +++ b/src/Models/HydrostaticFreeSurfaceModels/cache_hydrostatic_free_surface_tendencies.jl @@ -54,4 +54,3 @@ function cache_previous_tendencies!(model::HydrostaticFreeSurfaceModel) return nothing end - diff --git a/src/Models/HydrostaticFreeSurfaceModels/compute_hydrostatic_flux_bcs.jl b/src/Models/HydrostaticFreeSurfaceModels/compute_hydrostatic_flux_bcs.jl index ccef1df831e..b7b24d52668 100644 --- a/src/Models/HydrostaticFreeSurfaceModels/compute_hydrostatic_flux_bcs.jl +++ b/src/Models/HydrostaticFreeSurfaceModels/compute_hydrostatic_flux_bcs.jl @@ -20,10 +20,10 @@ function TimeSteppers.compute_flux_bc_tendencies!(model::HydrostaticFreeSurfaceM arch = architecture(grid) velocities = model.velocities tracers = model.tracers - + args = (model.clock, fields(model), model.closure, model.buoyancy) - + # Velocity fields for i in (:u, :v) @apply_regionally compute_flux_bcs!(Gⁿ[i], velocities[i], arch, args) diff --git a/src/Models/HydrostaticFreeSurfaceModels/vertical_vorticity.jl b/src/Models/HydrostaticFreeSurfaceModels/vertical_vorticity.jl index e23cf385995..ecdbbf801c7 100644 --- a/src/Models/HydrostaticFreeSurfaceModels/vertical_vorticity.jl +++ b/src/Models/HydrostaticFreeSurfaceModels/vertical_vorticity.jl @@ -14,4 +14,3 @@ function vertical_vorticity(model::HydrostaticFreeSurfaceModel) u, v, w = model.velocities return KernelFunctionOperation{Face, Face, Center}(ζ₃ᶠᶠᶜ, model.grid, u, v) end - diff --git a/src/Models/LagrangianParticleTracking/drogued_dynamics.jl b/src/Models/LagrangianParticleTracking/drogued_dynamics.jl index d39079ec1fc..bc2e5d6c929 100644 --- a/src/Models/LagrangianParticleTracking/drogued_dynamics.jl +++ b/src/Models/LagrangianParticleTracking/drogued_dynamics.jl @@ -1,8 +1,8 @@ """ DroguedParticleDynamics(depths) -`DroguedParticleDynamics` goes in the `dynamics` slot of `LagrangianParticles` -and modifies their behaviour to mimic the behaviour of bouys which are +`DroguedParticleDynamics` goes in the `dynamics` slot of `LagrangianParticles` +and modifies their behaviour to mimic the behaviour of bouys which are drogued at `depths`. The particles remain at the their `z` position so the "measurment depth can be set", and then are advected in the `x` and `y` directions according to the velocity field at `depths`. @@ -35,7 +35,7 @@ struct DroguedParticleDynamics{DD} depths :: DD end -Adapt.adapt_structure(to, dbd::DroguedParticleDynamics) = +Adapt.adapt_structure(to, dbd::DroguedParticleDynamics) = DroguedParticleDynamics(adapt(to, dbd.depths)) @inline (::DroguedParticleDynamics)(args...) = nothing diff --git a/src/Models/LagrangianParticleTracking/update_lagrangian_particle_properties.jl b/src/Models/LagrangianParticleTracking/update_lagrangian_particle_properties.jl index a4815f59377..148a4681746 100644 --- a/src/Models/LagrangianParticleTracking/update_lagrangian_particle_properties.jl +++ b/src/Models/LagrangianParticleTracking/update_lagrangian_particle_properties.jl @@ -34,4 +34,3 @@ function update_lagrangian_particle_properties!(particles, model, Δt) return nothing end - diff --git a/src/Models/NonhydrostaticModels/NonhydrostaticModels.jl b/src/Models/NonhydrostaticModels/NonhydrostaticModels.jl index 43f9cd3ddb2..7f932014ce5 100644 --- a/src/Models/NonhydrostaticModels/NonhydrostaticModels.jl +++ b/src/Models/NonhydrostaticModels/NonhydrostaticModels.jl @@ -51,7 +51,7 @@ nonhydrostatic_pressure_solver(arch, grid, ::Nothing) = ConjugateGradientPoisson const IBGWithFFT = ImmersedBoundaryGrid{<:Any, <:Any, <:Any, <:Any, <:GridWithFFTSolver} nonhydrostatic_pressure_solver(arch, ibg::IBGWithFFT, ::Nothing) = naive_solver_with_warning(arch, ibg, nothing) nonhydrostatic_pressure_solver(arch, ibg::IBGWithFFT, fs) = naive_solver_with_warning(arch, ibg, fs) - + function naive_solver_with_warning(arch, ibg, free_surface) msg = """The FFT-based pressure_solver for NonhydrostaticModels on ImmersedBoundaryGrid is approximate and will probably produce velocity fields that are divergent @@ -122,4 +122,3 @@ include("compute_nonhydrostatic_tendencies.jl") include("compute_nonhydrostatic_buffer_tendencies.jl") end # module - diff --git a/src/Models/NonhydrostaticModels/background_fields.jl b/src/Models/NonhydrostaticModels/background_fields.jl index fb065800ea8..9cde7cab630 100644 --- a/src/Models/NonhydrostaticModels/background_fields.jl +++ b/src/Models/NonhydrostaticModels/background_fields.jl @@ -21,7 +21,7 @@ function background_tracer_fields(bg, tracer_names, grid, clock) regularize_background_field(Center, Center, Center, getindex(bg, c), grid, clock) : ZeroField() for c in tracer_names) - + return NamedTuple{tracer_names}(tracer_fields) end @@ -37,7 +37,7 @@ struct BackgroundFields{Q, U, C} end end -Adapt.adapt_structure(to, bf::BackgroundFields{Q}) where Q = +Adapt.adapt_structure(to, bf::BackgroundFields{Q}) where Q = BackgroundFields{Q}(adapt(to, bf.velocities), adapt(to, bf.tracers)) const BackgroundFieldsWithClosureFluxes = BackgroundFields{true} @@ -111,7 +111,7 @@ function regularize_background_field(LX, LY, LZ, field::AbstractField, grid, clo if location(field) != (LX, LY, LZ) throw(ArgumentError("Cannot use field at $(location(field)) as a background field at $((LX, LY, LZ))")) end - + return field end diff --git a/src/Models/NonhydrostaticModels/compute_nonhydrostatic_buffer_tendencies.jl b/src/Models/NonhydrostaticModels/compute_nonhydrostatic_buffer_tendencies.jl index 7254a9e6d82..697d2cfac1e 100644 --- a/src/Models/NonhydrostaticModels/compute_nonhydrostatic_buffer_tendencies.jl +++ b/src/Models/NonhydrostaticModels/compute_nonhydrostatic_buffer_tendencies.jl @@ -81,4 +81,3 @@ function buffer_parameters(parameters, grid, arch) return Tuple(KernelParameters(parameters[i]...) for i in findall(include_side)) end - diff --git a/src/Models/NonhydrostaticModels/nonhydrostatic_tendency_kernel_functions.jl b/src/Models/NonhydrostaticModels/nonhydrostatic_tendency_kernel_functions.jl index 5d5ab18fb5a..a8220de003e 100644 --- a/src/Models/NonhydrostaticModels/nonhydrostatic_tendency_kernel_functions.jl +++ b/src/Models/NonhydrostaticModels/nonhydrostatic_tendency_kernel_functions.jl @@ -86,7 +86,7 @@ pressure anomaly. total_velocities = sum_of_velocities(velocities, background_fields.velocities) total_velocities = with_advective_forcing(forcing, total_velocities) - closure_velocities = assemble_closure_velocities(velocities, background_fields) + closure_velocities = assemble_closure_velocities(velocities, background_fields) closure_model_fields = merge(closure_velocities, tracers, auxiliary_fields) model_fields = merge(velocities, tracers, auxiliary_fields) @@ -148,7 +148,7 @@ pressure anomaly. total_velocities = sum_of_velocities(velocities, background_fields.velocities) total_velocities = with_advective_forcing(forcing, total_velocities) - closure_velocities = assemble_closure_velocities(velocities, background_fields) + closure_velocities = assemble_closure_velocities(velocities, background_fields) closure_model_fields = merge(closure_velocities, tracers, auxiliary_fields) model_fields = merge(velocities, tracers, auxiliary_fields) @@ -213,7 +213,7 @@ velocity components, tracer fields, and precalculated diffusivities where applic total_velocities = sum_of_velocities(velocities, background_fields.velocities) total_velocities = with_advective_forcing(forcing, total_velocities) - closure_velocities = assemble_closure_velocities(velocities, background_fields) + closure_velocities = assemble_closure_velocities(velocities, background_fields) closure_model_fields = merge(closure_velocities, tracers, auxiliary_fields) model_fields = merge(velocities, tracers, auxiliary_fields) @@ -296,4 +296,3 @@ velocity components, tracer fields, and precalculated diffusivities where applic + biogeochemical_transition(i, j, k, grid, biogeochemistry, val_tracer_name, clock, model_fields) + forcing(i, j, k, grid, clock, model_fields)) end - diff --git a/src/Models/NonhydrostaticModels/pressure_field.jl b/src/Models/NonhydrostaticModels/pressure_field.jl index 89d08ce898a..c7c0013c7d3 100644 --- a/src/Models/NonhydrostaticModels/pressure_field.jl +++ b/src/Models/NonhydrostaticModels/pressure_field.jl @@ -1,3 +1,2 @@ PressureField(model, data=model.pressures.pHY′.data; kw...) = Field(model.pressures.pHY′ + model.pressures.pNHS; data, kw...) - diff --git a/src/Models/NonhydrostaticModels/update_hydrostatic_pressure.jl b/src/Models/NonhydrostaticModels/update_hydrostatic_pressure.jl index 744425fd28b..be761ef64aa 100644 --- a/src/Models/NonhydrostaticModels/update_hydrostatic_pressure.jl +++ b/src/Models/NonhydrostaticModels/update_hydrostatic_pressure.jl @@ -49,4 +49,3 @@ update_hydrostatic_pressure!(::Nothing, arch, ::PCBIBG, args...; kw...) = nothin return KernelParameters(ii, jj) end - diff --git a/src/Models/ShallowWaterModels/compute_shallow_water_tendencies.jl b/src/Models/ShallowWaterModels/compute_shallow_water_tendencies.jl index d9f3c4138fc..38f1b2d635c 100644 --- a/src/Models/ShallowWaterModels/compute_shallow_water_tendencies.jl +++ b/src/Models/ShallowWaterModels/compute_shallow_water_tendencies.jl @@ -183,7 +183,7 @@ end """ Apply boundary conditions by adding flux divergences to the right-hand-side. """ function compute_flux_bc_tendencies!(model::ShallowWaterModel) - + Gⁿ = model.timestepper.Gⁿ arch = model.architecture clock = model.clock @@ -199,4 +199,3 @@ function compute_flux_bc_tendencies!(model::ShallowWaterModel) return nothing end - diff --git a/src/Models/ShallowWaterModels/store_shallow_water_tendencies.jl b/src/Models/ShallowWaterModels/store_shallow_water_tendencies.jl index 6add5b45632..f1d114d87eb 100644 --- a/src/Models/ShallowWaterModels/store_shallow_water_tendencies.jl +++ b/src/Models/ShallowWaterModels/store_shallow_water_tendencies.jl @@ -34,4 +34,3 @@ function cache_previous_tendencies!(model::ShallowWaterModel) return nothing end - diff --git a/src/Models/VarianceDissipationComputations/advective_dissipation.jl b/src/Models/VarianceDissipationComputations/advective_dissipation.jl index be9ad29d920..9480830c565 100644 --- a/src/Models/VarianceDissipationComputations/advective_dissipation.jl +++ b/src/Models/VarianceDissipationComputations/advective_dissipation.jl @@ -88,4 +88,4 @@ end Fⁿ.y[i, j, k] = _advective_tracer_flux_y(i, j, k, grid, advection, U.v, c) * σⁿ(i, j, k, grid, c, f, c) Fⁿ.z[i, j, k] = _advective_tracer_flux_z(i, j, k, grid, advection, U.w, c) * σⁿ(i, j, k, grid, c, c, f) end -end \ No newline at end of file +end diff --git a/src/Models/VarianceDissipationComputations/diffusive_dissipation.jl b/src/Models/VarianceDissipationComputations/diffusive_dissipation.jl index 2bcb46f5fd3..34fa3d0fe32 100644 --- a/src/Models/VarianceDissipationComputations/diffusive_dissipation.jl +++ b/src/Models/VarianceDissipationComputations/diffusive_dissipation.jl @@ -6,10 +6,10 @@ using Oceananigans.TurbulenceClosures: ExplicitTimeDiscretization C₁ = convert(eltype(grid), 1.5 + χ) C₂ = convert(eltype(grid), 0.5 + χ) - δˣc★ = δxᶠᶜᶜ(i, j, k, grid, c★, cⁿ⁺¹, cⁿ) + δˣc★ = δxᶠᶜᶜ(i, j, k, grid, c★, cⁿ⁺¹, cⁿ) δʸc★ = δyᶜᶠᶜ(i, j, k, grid, c★, cⁿ⁺¹, cⁿ) δᶻc★ = δzᶜᶜᶠ(i, j, k, grid, c★, cⁿ⁺¹, cⁿ) - + @inbounds begin vx₁ = C₁ * Vⁿ.x[i, j, k] / σⁿ(i, j, k, grid, f, c, c) vy₁ = C₁ * Vⁿ.y[i, j, k] / σⁿ(i, j, k, grid, c, f, c) @@ -28,58 +28,58 @@ end @kernel function _assemble_rk3_diffusive_dissipation!(K, grid, Vⁿ, cⁿ⁺¹, cⁿ) i, j, k = @index(Global, NTuple) - δˣc★ = δxᶠᶜᶜ(i, j, k, grid, c★, cⁿ⁺¹, cⁿ) + δˣc★ = δxᶠᶜᶜ(i, j, k, grid, c★, cⁿ⁺¹, cⁿ) δʸc★ = δyᶜᶠᶜ(i, j, k, grid, c★, cⁿ⁺¹, cⁿ) δᶻc★ = δzᶜᶜᶠ(i, j, k, grid, c★, cⁿ⁺¹, cⁿ) - + @inbounds begin vx₁ = Vⁿ.x[i, j, k] / σⁿ(i, j, k, grid, f, c, c) vy₁ = Vⁿ.y[i, j, k] / σⁿ(i, j, k, grid, c, f, c) vz₁ = Vⁿ.z[i, j, k] / σⁿ(i, j, k, grid, c, c, f) - K.x[i, j, k] = 2 * δˣc★ * vx₁ - K.y[i, j, k] = 2 * δʸc★ * vy₁ - K.z[i, j, k] = 2 * δᶻc★ * vz₁ + K.x[i, j, k] = 2 * δˣc★ * vx₁ + K.y[i, j, k] = 2 * δʸc★ * vy₁ + K.z[i, j, k] = 2 * δᶻc★ * vz₁ end end # QAB2 Implementation -@kernel function _cache_diffusive_fluxes!(Vⁿ, Vⁿ⁻¹, grid::AbstractGrid, clo, K, b, c, c_id, clk, fields) +@kernel function _cache_diffusive_fluxes!(Vⁿ, Vⁿ⁻¹, grid::AbstractGrid, clo, K, b, c, c_id, clk, fields) i, j, k = @index(Global, NTuple) - Vⁿ⁻¹.x[i, j, k] = Vⁿ.x[i, j, k] - Vⁿ⁻¹.y[i, j, k] = Vⁿ.y[i, j, k] - Vⁿ⁻¹.z[i, j, k] = Vⁿ.z[i, j, k] + Vⁿ⁻¹.x[i, j, k] = Vⁿ.x[i, j, k] + Vⁿ⁻¹.y[i, j, k] = Vⁿ.y[i, j, k] + Vⁿ⁻¹.z[i, j, k] = Vⁿ.z[i, j, k] Vⁿ.x[i, j, k] = zero(grid) Vⁿ.y[i, j, k] = zero(grid) Vⁿ.z[i, j, k] = zero(grid) - compute_diffusive_fluxes!(Vⁿ, i, j, k, grid, clo, K, b, c, c_id, clk, fields) + compute_diffusive_fluxes!(Vⁿ, i, j, k, grid, clo, K, b, c, c_id, clk, fields) end # RK3 Implementation for the last substep -@kernel function _cache_diffusive_fluxes!(Vⁿ, grid::AbstractGrid, clo, K, b, c, c_id, clk, fields) - i, j, k = @index(Global, NTuple) +@kernel function _cache_diffusive_fluxes!(Vⁿ, grid::AbstractGrid, clo, K, b, c, c_id, clk, fields) + i, j, k = @index(Global, NTuple) Vⁿ.x[i, j, k] = zero(grid) Vⁿ.y[i, j, k] = zero(grid) Vⁿ.z[i, j, k] = zero(grid) - compute_diffusive_fluxes!(Vⁿ, i, j, k, grid, clo, K, b, c, c_id, clk, fields) + compute_diffusive_fluxes!(Vⁿ, i, j, k, grid, clo, K, b, c, c_id, clk, fields) end # Deal with tuples of closures and diffusivities -@inline compute_diffusive_fluxes!(Vⁿ, i, j, k, grid, clo::Tuple{<:Any}, K, args...) = +@inline compute_diffusive_fluxes!(Vⁿ, i, j, k, grid, clo::Tuple{<:Any}, K, args...) = compute_diffusive_fluxes!(Vⁿ, i, j, k, grid, clo[1], K[1], args...) -@inline function compute_diffusive_fluxes!(Vⁿ, i, j, k, grid, clo::Tuple{<:Any, <:Any}, K, args...) +@inline function compute_diffusive_fluxes!(Vⁿ, i, j, k, grid, clo::Tuple{<:Any, <:Any}, K, args...) compute_diffusive_fluxes!(Vⁿ, i, j, k, grid, clo[1], K[1], args...) compute_diffusive_fluxes!(Vⁿ, i, j, k, grid, clo[2], K[2], args...) return nothing end -@inline function compute_diffusive_fluxes!(Vⁿ, i, j, k, grid, clo::Tuple, K, args...) +@inline function compute_diffusive_fluxes!(Vⁿ, i, j, k, grid, clo::Tuple, K, args...) compute_diffusive_fluxes!(Vⁿ, i, j, k, grid, clo[1], K[1], args...) compute_diffusive_fluxes!(Vⁿ, i, j, k, grid, clo[2], K[2], args...) compute_diffusive_fluxes!(Vⁿ, i, j, k, grid, clo[3:end], K[3:end], args...) diff --git a/src/Models/VarianceDissipationComputations/flatten_dissipation_fields.jl b/src/Models/VarianceDissipationComputations/flatten_dissipation_fields.jl index 8b8c55b675c..f40ce036d99 100644 --- a/src/Models/VarianceDissipationComputations/flatten_dissipation_fields.jl +++ b/src/Models/VarianceDissipationComputations/flatten_dissipation_fields.jl @@ -7,10 +7,10 @@ Flatten the dissipation fields of a `VarianceDissipation` object into a named tu - The dissipation associated with the closures in fields names `D-tracername-dir` - The squared gradients (necessary for computing an "effective diffusivity") in fields named `G-tracername-dir` """ -function flatten_dissipation_fields(t::VarianceDissipation) +function flatten_dissipation_fields(t::VarianceDissipation) A = t.advective_production D = t.diffusive_production - + tracer_name = t.tracer_name dirs = (:x, :y, :z) diff --git a/src/Models/VarianceDissipationComputations/update_fluxes.jl b/src/Models/VarianceDissipationComputations/update_fluxes.jl index f53d0edeb2a..c0e0eaa8f39 100644 --- a/src/Models/VarianceDissipationComputations/update_fluxes.jl +++ b/src/Models/VarianceDissipationComputations/update_fluxes.jl @@ -125,4 +125,4 @@ end @inbounds Uⁿ.u[i, j, k] = U.u[i, j, k] * Axᶠᶜᶜ(i, j, k, grid) @inbounds Uⁿ.v[i, j, k] = U.v[i, j, k] * Ayᶜᶠᶜ(i, j, k, grid) @inbounds Uⁿ.w[i, j, k] = U.w[i, j, k] * Azᶜᶜᶠ(i, j, k, grid) -end \ No newline at end of file +end diff --git a/src/Models/interleave_communication_and_computation.jl b/src/Models/interleave_communication_and_computation.jl index 53b11a85872..911ee304b86 100644 --- a/src/Models/interleave_communication_and_computation.jl +++ b/src/Models/interleave_communication_and_computation.jl @@ -65,4 +65,3 @@ function interior_tendency_kernel_parameters(arch::AsynchronousDistributed, grid return KernelParameters(sizes, offsets) end - diff --git a/src/MultiRegion/multi_region_models.jl b/src/MultiRegion/multi_region_models.jl index 4ccece6ebbb..15ada88f634 100644 --- a/src/MultiRegion/multi_region_models.jl +++ b/src/MultiRegion/multi_region_models.jl @@ -53,9 +53,9 @@ for T in Types end # TODO: For the moment, buoyancy gradients cannot be precomputed in MultiRegionModels -function BuoyancyForce(grid::MultiRegionGrids, formulation::AbstractBuoyancyFormulation; - gravity_unit_vector=NegativeZDirection(), - materialize_gradients=false) +function BuoyancyForce(grid::MultiRegionGrids, formulation::AbstractBuoyancyFormulation; + gravity_unit_vector=NegativeZDirection(), + materialize_gradients=false) gravity_unit_vector = validate_unit_vector(gravity_unit_vector) return BuoyancyForce(formulation, gravity_unit_vector, nothing) diff --git a/src/MultiRegion/multi_region_output_writers.jl b/src/MultiRegion/multi_region_output_writers.jl index 4cc08aa0a49..005e0cf1e1c 100644 --- a/src/MultiRegion/multi_region_output_writers.jl +++ b/src/MultiRegion/multi_region_output_writers.jl @@ -67,4 +67,3 @@ end function serializeproperty!(file, location, csg::ConformalCubedSphereGridOfSomeKind) file[location] = on_architecture(CPU(), csg) end - diff --git a/src/Operators/interpolation_operators.jl b/src/Operators/interpolation_operators.jl index b1f0e979c0c..deed2f885b9 100644 --- a/src/Operators/interpolation_operators.jl +++ b/src/Operators/interpolation_operators.jl @@ -153,4 +153,3 @@ end mask = active_nodes == 0 return ifelse(mask, zero(grid), ℑyzᵃᶜᶜ(i, j, k, grid, q, args...) / active_nodes) end - diff --git a/src/Operators/reciprocal_metric_operators.jl b/src/Operators/reciprocal_metric_operators.jl index b78c1916362..06e5fc0a6c9 100644 --- a/src/Operators/reciprocal_metric_operators.jl +++ b/src/Operators/reciprocal_metric_operators.jl @@ -8,8 +8,8 @@ for L1 in (:ᶜ, :ᶠ), L2 in (:ᶜ, :ᶠ, :ᵃ), L3 in (:ᶜ, :ᶠ, :ᵃ) metric_x = Symbol(func, :x, L1, L2, L3) metric_y = Symbol(func, :y, L2, L1, L3) metric_z = Symbol(func, :z, L2, L3, L1) - - @eval begin + + @eval begin @inline $rcp_x(i, j, k, grid) = 1 / $metric_x(i, j, k, grid) @inline $rcp_y(i, j, k, grid) = 1 / $metric_y(i, j, k, grid) @inline $rcp_z(i, j, k, grid) = 1 / $metric_z(i, j, k, grid) @@ -24,8 +24,7 @@ for L1 in (:ᶜ, :ᶠ), L2 in (:ᶜ, :ᶠ), L3 in (:ᶜ, :ᶠ) volume = Symbol(:V, L1, L2, L3) @eval begin @inline $rcp_volume(i, j, k, grid) = 1 / $volume(i, j, k, grid) - + export $rcp_volume end end - diff --git a/src/Operators/spacings_and_areas_and_volumes.jl b/src/Operators/spacings_and_areas_and_volumes.jl index 00cbfe7d5a0..dec6fbb9cda 100644 --- a/src/Operators/spacings_and_areas_and_volumes.jl +++ b/src/Operators/spacings_and_areas_and_volumes.jl @@ -217,8 +217,8 @@ end for sym in (:Δxᶠᶜᵃ, :Δxᶜᶠᵃ, :Δxᶠᶠᵃ, :Δxᶜᶜᵃ) @eval @inline function $sym(i::AbstractArray, j::AbstractArray, k::AbstractArray, grid::LLGX) - x = $sym(1, j, 1, grid) - Base.Broadcast.materialize(Base.Broadcast.Broadcasted(Base.Broadcast.BroadcastStyle(typeof(x)), Base.identity, (Base.Broadcast.Extruded(x, (false, true,false), (1,1,1)),), (Base.axes(i, 1), Base.axes(j, 2), Base.axes(k, 3)))) + x = $sym(1, j, 1, grid) + Base.Broadcast.materialize(Base.Broadcast.Broadcasted(Base.Broadcast.BroadcastStyle(typeof(x)), Base.identity, (Base.Broadcast.Extruded(x, (false, true,false), (1,1,1)),), (Base.axes(i, 1), Base.axes(j, 2), Base.axes(k, 3)))) end end @@ -382,12 +382,12 @@ end ##### We also use the function "volume" rather than `V`. ##### -function location_from_superscript(val::Symbol) - if val == :ᶜ - return :Center - elseif val == :ᶠ - return :Face - else +function location_from_superscript(val::Symbol) + if val == :ᶜ + return :Center + elseif val == :ᶠ + return :Face + else return :Nothing end end @@ -401,10 +401,10 @@ for ℓ1 in (:ᶜ, :ᶠ), ℓ2 in (:ᶜ, :ᶠ, :ᵃ), ℓ3 in (:ᶜ, :ᶠ, :ᵃ) L3 = location_from_superscript(ℓ3) spacing_x = Symbol(:Δx, ℓ1, ℓ2, ℓ3) - spacing_λ = Symbol(:Δλ, ℓ1, ℓ2, ℓ3) + spacing_λ = Symbol(:Δλ, ℓ1, ℓ2, ℓ3) spacing_y = Symbol(:Δy, ℓ2, ℓ1, ℓ3) - spacing_φ = Symbol(:Δφ, ℓ2, ℓ1, ℓ3) - spacing_z = Symbol(:Δz, ℓ2, ℓ3, ℓ1) + spacing_φ = Symbol(:Δφ, ℓ2, ℓ1, ℓ3) + spacing_z = Symbol(:Δz, ℓ2, ℓ3, ℓ1) spacing_r = Symbol(:Δr, ℓ2, ℓ3, ℓ1) @eval begin @@ -428,7 +428,7 @@ for ℓ1 in (:ᶜ, :ᶠ), ℓ2 in (:ᶜ, :ᶠ), ℓ3 in (:ᶜ, :ᶠ, :ᵃ) area_x = Symbol(:Ax, ℓ3, ℓ1, ℓ2) area_y = Symbol(:Ay, ℓ1, ℓ3, ℓ2) area_z = Symbol(:Az, ℓ1, ℓ2, ℓ3) - + @eval begin @eval Ax(i, j, k, grid, ::$L3, ::$L1, ::$L2) = $area_x(i, j, k, grid) @eval Ay(i, j, k, grid, ::$L1, ::$L3, ::$L2) = $area_y(i, j, k, grid) diff --git a/src/OrthogonalSphericalShellGrids/distributed_zipper.jl b/src/OrthogonalSphericalShellGrids/distributed_zipper.jl index c20b928fdff..ebd5fa32246 100644 --- a/src/OrthogonalSphericalShellGrids/distributed_zipper.jl +++ b/src/OrthogonalSphericalShellGrids/distributed_zipper.jl @@ -59,7 +59,7 @@ end end function fill_halo_regions!(c::OffsetArray, bcs, indices, loc, grid::DistributedTripolarGridOfSomeKind, buffers, args...; kwargs...) - + arch = architecture(grid) kernels!, ordered_bcs = get_boundary_kernels(bcs, c, grid, loc, indices) @@ -83,7 +83,7 @@ function fill_halo_regions!(c::OffsetArray, bcs, indices, loc, grid::Distributed north_bc = bcs.north switch_north_halos!(c, north_bc, grid, loc) end - + return nothing end diff --git a/src/OutputReaders/OutputReaders.jl b/src/OutputReaders/OutputReaders.jl index 6df31e1d88e..04b3dac221d 100644 --- a/src/OutputReaders/OutputReaders.jl +++ b/src/OutputReaders/OutputReaders.jl @@ -41,4 +41,3 @@ include("combining_field_time_series.jl") include("field_dataset.jl") end # module - diff --git a/src/OutputReaders/show_field_time_series.jl b/src/OutputReaders/show_field_time_series.jl index 8251a0ec976..4c7e3d3ce2d 100644 --- a/src/OutputReaders/show_field_time_series.jl +++ b/src/OutputReaders/show_field_time_series.jl @@ -58,4 +58,3 @@ field_time_series_suffix(fts::OnDiskFTS) = string("├── backend: ", summary(fts.backend), '\n', "├── path: ", fts.path, '\n', "└── name: ", fts.name) - diff --git a/src/Simulations/time_step_wizard.jl b/src/Simulations/time_step_wizard.jl index 8b284cf780a..7afe0f21c42 100644 --- a/src/Simulations/time_step_wizard.jl +++ b/src/Simulations/time_step_wizard.jl @@ -127,4 +127,3 @@ function conjure_time_step_wizard!(simulation, schedule=IterationInterval(10); w simulation.callbacks[:time_step_wizard] = Callback(wizard, schedule) return nothing end - diff --git a/src/Solvers/poisson_eigenvalues.jl b/src/Solvers/poisson_eigenvalues.jl index 8a0dba0f87f..b36d6ddb941 100644 --- a/src/Solvers/poisson_eigenvalues.jl +++ b/src/Solvers/poisson_eigenvalues.jl @@ -29,4 +29,3 @@ Return N-element array of `0.0` reshaped to three-dimensions. This is also the first `poisson_eigenvalue` for `Bounded` and `Periodic` directions. """ poisson_eigenvalues(N, L, dim, ::Flat) = reshape(zeros(N), reshaped_size(N, dim)...) - diff --git a/src/TurbulenceClosures/abstract_scalar_diffusivity_closure.jl b/src/TurbulenceClosures/abstract_scalar_diffusivity_closure.jl index 100b6b81433..22d32cb38ee 100644 --- a/src/TurbulenceClosures/abstract_scalar_diffusivity_closure.jl +++ b/src/TurbulenceClosures/abstract_scalar_diffusivity_closure.jl @@ -112,7 +112,7 @@ const AVD = AbstractScalarDiffusivity{<:Any, <:VerticalFormulation} @inline κᶠᶜᶠ(i, j, k, grid, closure::ASD, K, id, clk, fields) = κᶠᶜᶠ(i, j, k, grid, diffusivity_location(closure), diffusivity(closure, K, id), clk, fields) @inline κᶜᶠᶠ(i, j, k, grid, closure::ASD, K, id, clk, fields) = κᶜᶠᶠ(i, j, k, grid, diffusivity_location(closure), diffusivity(closure, K, id), clk, fields) -# Vertical and horizontal diffusivity +# Vertical and horizontal diffusivity # Viscosities without explicit passing of `id` @inline νzᶜᶜᶜ(i, j, k, grid, closure::ASD, K, clk, fields) = νᶜᶜᶜ(i, j, k, grid, closure, K, clk, fields) @@ -285,7 +285,7 @@ end zero(grid)) end -@inline function diffusive_flux_z(i, j, k, grid::VerticallyBoundedGrid, ::VITD, closure::AIDorAVD, K, id, c, clk, fields, b) +@inline function diffusive_flux_z(i, j, k, grid::VerticallyBoundedGrid, ::VITD, closure::AIDorAVD, K, id, c, clk, fields, b) return ifelse((k == 1) | (k == grid.Nz+1), diffusive_flux_z(i, j, k, grid, ExplicitTimeDiscretization(), closure, K, id, c, clk, fields, b), zero(grid)) diff --git a/src/TurbulenceClosures/closure_fields.jl b/src/TurbulenceClosures/closure_fields.jl index ee7bff699f3..ce39684a5a1 100644 --- a/src/TurbulenceClosures/closure_fields.jl +++ b/src/TurbulenceClosures/closure_fields.jl @@ -25,4 +25,3 @@ build_closure_fields(grid, clock, tracer_names, bcs, closure) = nothing build_closure_fields(grid, clock, tracer_names, bcs, closure_tuple::Tuple) = Tuple(build_closure_fields(grid, clock, tracer_names, bcs, closure) for closure in closure_tuple) - diff --git a/src/TurbulenceClosures/turbulence_closure_implementations/Smagorinskys/lilly_coefficient.jl b/src/TurbulenceClosures/turbulence_closure_implementations/Smagorinskys/lilly_coefficient.jl index d183db66675..555c199c68d 100644 --- a/src/TurbulenceClosures/turbulence_closure_implementations/Smagorinskys/lilly_coefficient.jl +++ b/src/TurbulenceClosures/turbulence_closure_implementations/Smagorinskys/lilly_coefficient.jl @@ -145,5 +145,3 @@ Base.summary(dc::LillyCoefficient) = string("LillyCoefficient(smagorinsky = $(dc Base.show(io::IO, dc::LillyCoefficient) = print(io, "LillyCoefficient with\n", "├── Smagorinsky coefficient = ", dc.smagorinsky, "\n", "└── reduction_factor = ", dc.reduction_factor) - - diff --git a/src/TurbulenceClosures/turbulence_closure_implementations/TKEBasedVerticalDiffusivities/time_step_catke_equation.jl b/src/TurbulenceClosures/turbulence_closure_implementations/TKEBasedVerticalDiffusivities/time_step_catke_equation.jl index 056c5d4f6b4..4503efaf937 100644 --- a/src/TurbulenceClosures/turbulence_closure_implementations/TKEBasedVerticalDiffusivities/time_step_catke_equation.jl +++ b/src/TurbulenceClosures/turbulence_closure_implementations/TKEBasedVerticalDiffusivities/time_step_catke_equation.jl @@ -57,7 +57,7 @@ function time_step_catke_equation!(model, ::QuasiAdamsBashforth2TimeStepper) compute_TKE_diffusivity!, κe, grid, closure, model.velocities, tracers, buoyancy, closure_fields) - + # ... and step forward. launch!(arch, grid, :xyz, _ab2_substep_turbulent_kinetic_energy!, @@ -74,8 +74,8 @@ function time_step_catke_equation!(model, ::QuasiAdamsBashforth2TimeStepper) implicit_step!(e, implicit_solver, closure, closure_fields, Val(tracer_index), - model.clock, - fields(model), + model.clock, + fields(model), Δτ) end @@ -84,8 +84,8 @@ end @inline rk3_coeffs(ts, stage) = stage == 1 ? (one(ts.γ²), zero(ts.γ²)) : stage == 2 ? (ts.γ², ts.ζ²) : - (ts.γ³, ts.ζ³) - + (ts.γ³, ts.ζ³) + function time_step_catke_equation!(model, ::SplitRungeKutta3TimeStepper) # TODO: properly handle closure tuples @@ -121,7 +121,7 @@ function time_step_catke_equation!(model, ::SplitRungeKutta3TimeStepper) compute_TKE_diffusivity!, κe, grid, closure, model.velocities, tracers, buoyancy, closure_fields) - + # ... and step forward. launch!(arch, grid, :xyz, _euler_step_turbulent_kinetic_energy!, @@ -132,10 +132,10 @@ function time_step_catke_equation!(model, ::SplitRungeKutta3TimeStepper) implicit_step!(e, implicit_solver, closure, closure_fields, Val(tracer_index), - model.clock, - fields(model), + model.clock, + fields(model), Δt) - + return nothing end @@ -242,7 +242,7 @@ end Δτ = convert(FT, Δτ) e = tracers.e - # See below. + # See below. α = convert(FT, 1.5) + χ β = convert(FT, 0.5) + χ diff --git a/src/TurbulenceClosures/turbulence_closure_implementations/TKEBasedVerticalDiffusivities/tke_dissipation_stability_functions.jl b/src/TurbulenceClosures/turbulence_closure_implementations/TKEBasedVerticalDiffusivities/tke_dissipation_stability_functions.jl index df281373c3b..bf9c8512c01 100644 --- a/src/TurbulenceClosures/turbulence_closure_implementations/TKEBasedVerticalDiffusivities/tke_dissipation_stability_functions.jl +++ b/src/TurbulenceClosures/turbulence_closure_implementations/TKEBasedVerticalDiffusivities/tke_dissipation_stability_functions.jl @@ -291,4 +291,3 @@ end return num / den end - diff --git a/src/TurbulenceClosures/turbulence_closure_implementations/advective_skew_diffusion.jl b/src/TurbulenceClosures/turbulence_closure_implementations/advective_skew_diffusion.jl index abccbe2f00e..6c696818779 100644 --- a/src/TurbulenceClosures/turbulence_closure_implementations/advective_skew_diffusion.jl +++ b/src/TurbulenceClosures/turbulence_closure_implementations/advective_skew_diffusion.jl @@ -27,7 +27,7 @@ end tapering_factor(Sx, Sy, slope_limiter) Return the tapering factor `min(1, Sₘₐₓ² / S²)`, where `S² = Sx² + Sy²` -that multiplies all components of the isopycnal slope tensor. +that multiplies all components of the isopycnal slope tensor. References ========== @@ -36,27 +36,27 @@ R. Gerdes, C. Koberle, and J. Willebrand. (1991), "The influence of numerical ad """ @inline function tapering_factor(Sx, Sy, slope_limiter::FluxTapering) S² = Sx^2 + Sy^2 - Sₘ² = slope_limiter.max_slope^2 + Sₘ² = slope_limiter.max_slope^2 return min(one(Sx), Sₘ² / S²) end # Slope in x-direction at F, C, F locations, zeroed out on peripheries -@inline function Sxᶠᶜᶠ(i, j, k, grid, b, C) - bx = ℑzᵃᵃᶠ(i, j, k, grid, ∂x_b, b, C) +@inline function Sxᶠᶜᶠ(i, j, k, grid, b, C) + bx = ℑzᵃᵃᶠ(i, j, k, grid, ∂x_b, b, C) bz = ℑxᶠᵃᵃ(i, j, k, grid, ∂z_b, b, C) - + Sx = ifelse(bz == 0, zero(grid), - bx / bz) # Impose a boundary condition on immersed peripheries inactive = peripheral_node(i, j, k, grid, Face(), Center(), Face()) Sx = ifelse(inactive, zero(grid), Sx) - + return Sx end # Slope in y-direction at F, C, F locations, zeroed out on peripheries -@inline function Syᶜᶠᶠ(i, j, k, grid, b, C) - by = ℑzᵃᵃᶠ(i, j, k, grid, ∂y_b, b, C) +@inline function Syᶜᶠᶠ(i, j, k, grid, b, C) + by = ℑzᵃᵃᶠ(i, j, k, grid, ∂y_b, b, C) bz = ℑyᵃᶠᵃ(i, j, k, grid, ∂z_b, b, C) Sy = ifelse(bz == 0, zero(grid), - by / bz) @@ -64,20 +64,20 @@ end # Impose a boundary condition on immersed peripheries inactive = peripheral_node(i, j, k, grid, Center(), Face(), Face()) Sy = ifelse(inactive, zero(grid), Sy) - + return Sy end # tapered slope in x-direction at F, C, F locations @inline function ϵSxᶠᶜᶠ(i, j, k, grid, slope_limiter, b, C) - Sx = Sxᶠᶜᶠ(i, j, k, grid, b, C) + Sx = Sxᶠᶜᶠ(i, j, k, grid, b, C) ϵ = tapering_factor(Sx, zero(grid), slope_limiter) return ϵ * Sx end # tapered slope in y-direction at F, C, F locations @inline function ϵSyᶜᶠᶠ(i, j, k, grid, slope_limiter, b, C) - Sy = Syᶜᶠᶠ(i, j, k, grid, b, C) + Sy = Syᶜᶠᶠ(i, j, k, grid, b, C) ϵ = tapering_factor(zero(grid), Sy, slope_limiter) return ϵ * Sy end @@ -96,9 +96,9 @@ end uₑ[i, j, k] = - δzᵃᵃᶜ(i, j, k, grid, κ_ϵSxᶠᶜᶠ, clock, slope_limiter, κ, buoyancy, fields) * Δz⁻¹ᶠᶜᶜ(i, j, k, grid) vₑ[i, j, k] = - δzᵃᵃᶜ(i, j, k, grid, κ_ϵSyᶜᶠᶠ, clock, slope_limiter, κ, buoyancy, fields) * Δz⁻¹ᶜᶠᶜ(i, j, k, grid) - wˣ = δxᶜᵃᵃ(i, j, k, grid, Δy_qᶠᶜᶠ, κ_ϵSxᶠᶜᶠ, clock, slope_limiter, κ, buoyancy, fields) - wʸ = δyᵃᶜᵃ(i, j, k, grid, Δx_qᶜᶠᶠ, κ_ϵSyᶜᶠᶠ, clock, slope_limiter, κ, buoyancy, fields) - + wˣ = δxᶜᵃᵃ(i, j, k, grid, Δy_qᶠᶜᶠ, κ_ϵSxᶠᶜᶠ, clock, slope_limiter, κ, buoyancy, fields) + wʸ = δyᵃᶜᵃ(i, j, k, grid, Δx_qᶜᶠᶠ, κ_ϵSyᶜᶠᶠ, clock, slope_limiter, κ, buoyancy, fields) + wₑ[i, j, k] = (wˣ + wʸ) * Az⁻¹ᶜᶜᶠ(i, j, k, grid) end end diff --git a/src/TurbulenceClosures/turbulence_closure_utils.jl b/src/TurbulenceClosures/turbulence_closure_utils.jl index d52d480406b..4fe91c75c43 100644 --- a/src/TurbulenceClosures/turbulence_closure_utils.jl +++ b/src/TurbulenceClosures/turbulence_closure_utils.jl @@ -34,4 +34,3 @@ end return NamedTuple{κ_names}(κ_values) end - diff --git a/src/TurbulenceClosures/vertically_implicit_diffusion_solver.jl b/src/TurbulenceClosures/vertically_implicit_diffusion_solver.jl index da527def12b..93009fdd49a 100644 --- a/src/TurbulenceClosures/vertically_implicit_diffusion_solver.jl +++ b/src/TurbulenceClosures/vertically_implicit_diffusion_solver.jl @@ -115,13 +115,13 @@ end # Fallback for single closure. These coefficients are extended for tupled closures in `closure_tuples.jl` -@inline _implicit_linear_coefficient(i, j, k, grid, closure, K, id, ℓx, ℓy, ℓz, Δt, clock, fields) = +@inline _implicit_linear_coefficient(i, j, k, grid, closure, K, id, ℓx, ℓy, ℓz, Δt, clock, fields) = implicit_linear_coefficient(i, j, k, grid, closure, K, id, ℓx, ℓy, ℓz, Δt, clock, fields) -@inline _ivd_upper_diagonal(i, j, k, grid, closure, K, id, ℓx, ℓy, ℓz, Δt, clock, fields) = +@inline _ivd_upper_diagonal(i, j, k, grid, closure, K, id, ℓx, ℓy, ℓz, Δt, clock, fields) = ivd_upper_diagonal(i, j, k, grid, closure, K, id, ℓx, ℓy, ℓz, Δt, clock, fields) -@inline _ivd_lower_diagonal(i, j, k, grid, closure, K, id, ℓx, ℓy, ℓz, Δt, clock, fields) = +@inline _ivd_lower_diagonal(i, j, k, grid, closure, K, id, ℓx, ℓy, ℓz, Δt, clock, fields) = ivd_lower_diagonal(i, j, k, grid, closure, K, id, ℓx, ℓy, ℓz, Δt, clock, fields) ##### @@ -166,13 +166,13 @@ end # Extend `get_coefficient` to retrieve `ivd_diagonal`, `_ivd_lower_diagonal` and `_ivd_upper_diagonal`. # Note that we use the "periphery-aware" upper and lower diagonals -@inline get_coefficient(i, j, k, grid, ::VerticallyImplicitDiffusionLowerDiagonal, p, ::ZDirection, clo, K, id, ℓx, ℓy, ℓz, Δt, clk, fields) = +@inline get_coefficient(i, j, k, grid, ::VerticallyImplicitDiffusionLowerDiagonal, p, ::ZDirection, clo, K, id, ℓx, ℓy, ℓz, Δt, clk, fields) = _ivd_lower_diagonal(i, j, k, grid, clo, K, id, ℓx, ℓy, ℓz, Δt, clk, fields) -@inline get_coefficient(i, j, k, grid, ::VerticallyImplicitDiffusionUpperDiagonal, p, ::ZDirection, clo, K, id, ℓx, ℓy, ℓz, Δt, clk, fields) = +@inline get_coefficient(i, j, k, grid, ::VerticallyImplicitDiffusionUpperDiagonal, p, ::ZDirection, clo, K, id, ℓx, ℓy, ℓz, Δt, clk, fields) = _ivd_upper_diagonal(i, j, k, grid, clo, K, id, ℓx, ℓy, ℓz, Δt, clk, fields) -@inline get_coefficient(i, j, k, grid, ::VerticallyImplicitDiffusionDiagonal, p, ::ZDirection, clo, K, id, ℓx, ℓy, ℓz, Δt, clk, fields) = +@inline get_coefficient(i, j, k, grid, ::VerticallyImplicitDiffusionDiagonal, p, ::ZDirection, clo, K, id, ℓx, ℓy, ℓz, Δt, clk, fields) = ivd_diagonal(i, j, k, grid, clo, K, id, ℓx, ℓy, ℓz, Δt, clk, fields) ##### @@ -196,9 +196,9 @@ function implicit_step!(field::Field, closure_fields, tracer_index, clock, - fields, + fields, Δt) - + # Filter explicit closures for closure tuples if closure isa Tuple closure_tuple = closure diff --git a/src/Utils/prettysummary.jl b/src/Utils/prettysummary.jl index 06de620d699..f12782243e8 100644 --- a/src/Utils/prettysummary.jl +++ b/src/Utils/prettysummary.jl @@ -43,4 +43,3 @@ function prettykeys(t) length(names) == 1 && return string(first(names)) return string("(", (string(n, ", ") for n in names[1:end-1])..., last(names), ")") end - diff --git a/src/Utils/sum_of_arrays.jl b/src/Utils/sum_of_arrays.jl index 22970811517..c56e8ba1a93 100644 --- a/src/Utils/sum_of_arrays.jl +++ b/src/Utils/sum_of_arrays.jl @@ -56,4 +56,4 @@ const NT = NamedTuple @inline sum_of_velocities(U1::NT, ::Nothing, ::Nothing) = U1 @inline sum_of_velocities(::Nothing, U2::NT, ::Nothing) = U2 -@inline sum_of_velocities(::Nothing, ::Nothing, U3::NT) = U3 \ No newline at end of file +@inline sum_of_velocities(::Nothing, ::Nothing, U3::NT) = U3 diff --git a/src/Utils/tabulated_function.jl b/src/Utils/tabulated_function.jl index 1c2f3d77f4e..d656be0870c 100644 --- a/src/Utils/tabulated_function.jl +++ b/src/Utils/tabulated_function.jl @@ -172,4 +172,3 @@ function Base.show(io::IO, f::TabulatedFunction) print(io, " of ", f.func) end end - diff --git a/src/Utils/user_function_arguments.jl b/src/Utils/user_function_arguments.jl index e17045a5494..f108944f838 100644 --- a/src/Utils/user_function_arguments.jl +++ b/src/Utils/user_function_arguments.jl @@ -37,4 +37,3 @@ end return tuple(field_args..., parameters) end - diff --git a/test/distributed_tests_utils.jl b/test/distributed_tests_utils.jl index d04a05af55e..ca49c5d334d 100644 --- a/test/distributed_tests_utils.jl +++ b/test/distributed_tests_utils.jl @@ -165,4 +165,3 @@ function run_distributed_simulation(grid) return model end - diff --git a/test/reactant_test_utils.jl b/test/reactant_test_utils.jl index 99c3d8a5205..a6feda311f8 100644 --- a/test/reactant_test_utils.jl +++ b/test/reactant_test_utils.jl @@ -87,7 +87,7 @@ function test_reactant_model_correctness(GridType, ModelType, grid_kw, model_kw; Oceananigans.TimeSteppers.update_state!(model) Oceananigans.Models.HydrostaticFreeSurfaceModels.compute_hydrostatic_momentum_tendencies!(model, model.velocities, :xyz) @jit Oceananigans.TimeSteppers.update_state!(r_model) - + mod = @code_hlo optimize=:before_jit Oceananigans.Models.HydrostaticFreeSurfaceModels.compute_hydrostatic_momentum_tendencies!(r_model, r_model.velocities, :xyz) @jit Oceananigans.Models.HydrostaticFreeSurfaceModels.compute_hydrostatic_momentum_tendencies!(r_model, r_model.velocities, :xyz) @@ -196,4 +196,3 @@ end i, j, k = @index(Global, NTuple) @inbounds f[i, j, k] += 1 end - diff --git a/test/test_background_flux_divergence.jl b/test/test_background_flux_divergence.jl index 361aac0909a..9d274fb52af 100644 --- a/test/test_background_flux_divergence.jl +++ b/test/test_background_flux_divergence.jl @@ -2,20 +2,20 @@ include("dependencies_for_runtests.jl") """ function run_with_background_fields(arch; with_background=true) - + Run a model with or without background fields and compare the two. """ function run_with_background_fields(arch; with_background=true) grid = RectilinearGrid(arch, size=10, z=(0, 1), topology=(Flat, Flat, Bounded)) # Setup model with or without background fields if with_background - background_fields = Oceananigans.BackgroundFields(; + background_fields = Oceananigans.BackgroundFields(; background_closure_fluxes=true, b=B̄_field) # we want no flux bottom boundary (∂B∂z = 0) and infinite ocean at the top boundary B_bcs = FieldBoundaryConditions( bottom = GradientBoundaryCondition(-N^2), # ∂B∂z = 0 → ∂b∂z = -∂B∂z = -N² - top = GradientBoundaryCondition(0.) # ∂B∂z = ∂B̄∂z+∂b∂z = N² → ∂b∂z = 0 - ); + top = GradientBoundaryCondition(0.) # ∂B∂z = ∂B̄∂z+∂b∂z = N² → ∂b∂z = 0 + ); model = NonhydrostaticModel(grid; background_fields, tracers = :b, buoyancy=BuoyancyTracer(), boundary_conditions=(; b = B_bcs)) b = model.tracers.b @@ -38,7 +38,7 @@ function run_with_background_fields(arch; with_background=true) # Run for a few iterations simulation = Simulation(model, Δt=0.1, stop_iteration=5) run!(simulation) - + return B end @@ -51,13 +51,13 @@ test_archs = has_cuda() ? [CPU(), GPU()] : [CPU()] @testset "Background Fields Tests" begin for arch in test_archs - + # Test model runs with background fields @test run_with_background_fields(arch, with_background=true) !== nothing - + # Test model runs without background fields @test run_with_background_fields(arch, with_background=false) !== nothing - + # Test that background fields affect the solution b_with = run_with_background_fields(arch, with_background=true) b_without = run_with_background_fields(arch, with_background=false) @@ -65,4 +65,4 @@ test_archs = has_cuda() ? [CPU(), GPU()] : [CPU()] # to pass the test, both total buoyancy should be the same even the method is not the same @test all(isapprox.(b_with, b_without, rtol=1e-10)) end -end \ No newline at end of file +end diff --git a/test/test_batched_tridiagonal_solver.jl b/test/test_batched_tridiagonal_solver.jl index 61d7e617e20..2699500e8ef 100644 --- a/test/test_batched_tridiagonal_solver.jl +++ b/test/test_batched_tridiagonal_solver.jl @@ -112,4 +112,3 @@ end end end end - diff --git a/test/test_computed_field.jl b/test/test_computed_field.jl index 5291e7245f3..4171a09d100 100644 --- a/test/test_computed_field.jl +++ b/test/test_computed_field.jl @@ -624,4 +624,3 @@ for arch in archs end end end - diff --git a/test/test_distributed_architectures.jl b/test/test_distributed_architectures.jl index 08303a71a42..601591e16a1 100644 --- a/test/test_distributed_architectures.jl +++ b/test/test_distributed_architectures.jl @@ -80,4 +80,3 @@ end end end =# - diff --git a/test/test_distributed_transpose.jl b/test/test_distributed_transpose.jl index 49b92e707eb..8f9167eba2a 100644 --- a/test/test_distributed_transpose.jl +++ b/test/test_distributed_transpose.jl @@ -52,4 +52,3 @@ end @test test_transpose((16, 44, 8), (2, 2, 1), topology, child_arch) end end - diff --git a/test/test_enzyme.jl b/test/test_enzyme.jl index d1f30ffca39..a7855824d17 100644 --- a/test/test_enzyme.jl +++ b/test/test_enzyme.jl @@ -378,4 +378,3 @@ end rel_error = abs(dedν[1][1] - ΔeΔν) / abs(ΔeΔν) @test rel_error < tol end - diff --git a/test/test_field.jl b/test/test_field.jl index 579cd6a837f..7fa4973f368 100644 --- a/test/test_field.jl +++ b/test/test_field.jl @@ -1,810 +1,810 @@ -include("dependencies_for_runtests.jl") - -using Statistics - -using Oceananigans.Fields: CenterField, ReducedField, has_velocities -using Oceananigans.Fields: VelocityFields, TracerFields, interpolate, interpolate! -using Oceananigans.Fields: reduced_location -using Oceananigans.Fields: FractionalIndices, interpolator, instantiate -using Oceananigans.Fields: convert_to_0_360, convert_to_λ₀_λ₀_plus360 -using Oceananigans.Grids: ξnode, ηnode, rnode -using Oceananigans.Grids: total_length -using Oceananigans.Grids: λnode -using Oceananigans.Grids: RectilinearGrid -using Oceananigans.ImmersedBoundaries: ImmersedBoundaryGrid, GridFittedBottom - -using Random -using GPUArraysCore: @allowscalar - -""" - correct_field_size(grid, FieldType, Tx, Ty, Tz) - -Test that the field initialized by the FieldType constructor on `grid` -has size `(Tx, Ty, Tz)`. -""" -correct_field_size(grid, loc, Tx, Ty, Tz) = size(parent(Field(instantiate(loc), grid))) == (Tx, Ty, Tz) - -function run_similar_field_tests(f) - g = similar(f) - @test typeof(f) == typeof(g) - @test f.grid == g.grid - @test location(f) === location(g) - @test !(f.data === g.data) - return nothing -end - -""" - correct_field_value_was_set(N, L, ftf, val) - -Test that the field initialized by the field type function `ftf` on the grid g -can be correctly filled with the value `val` using the `set!(f::AbstractField, v)` -function. -""" -function correct_field_value_was_set(grid, FieldType, val::Number) - arch = architecture(grid) - f = FieldType(grid) - set!(f, val) - return all(interior(f) .≈ val * on_architecture(arch, ones(size(f)))) -end - -function run_field_reduction_tests(grid) - u = XFaceField(grid) - v = YFaceField(grid) - w = ZFaceField(grid) - c = CenterField(grid) - η = Field{Center, Center, Nothing}(grid) - - f(x, y, z) = 1 + exp(x) * sin(y) * tanh(z) - - ϕs = [u, v, w, c] - [set!(ϕ, f) for ϕ in ϕs] - - z_top = znodes(grid, Face())[end] - set!(η, (x, y) -> f(x, y, z_top)) - push!(ϕs, η) - ϕs = Tuple(ϕs) - - u_vals = f.(nodes(u, reshape=true)...) - v_vals = f.(nodes(v, reshape=true)...) - w_vals = f.(nodes(w, reshape=true)...) - c_vals = f.(nodes(c, reshape=true)...) - η_vals = f.(nodes(η, reshape=true)...) - - # Convert to CuArray if needed. - arch = architecture(grid) - u_vals = on_architecture(arch, u_vals) - v_vals = on_architecture(arch, v_vals) - w_vals = on_architecture(arch, w_vals) - c_vals = on_architecture(arch, c_vals) - η_vals = on_architecture(arch, η_vals) - - ϕs_vals = (u_vals, v_vals, w_vals, c_vals, η_vals) - - dims_to_test = (1, 2, 3, (1, 2), (1, 3), (2, 3), (1, 2, 3)) - - for (ϕ, ϕ_vals) in zip(ϕs, ϕs_vals) - ε = eps(eltype(grid)) * 10 * maximum(maximum.(ϕs_vals)) - @info " Testing field reductions with tolerance $ε..." - - @test @allowscalar all(isapprox.(ϕ, ϕ_vals, atol=ε)) # if this isn't true, reduction tests can't pass - - # Important to make sure no scalar operations occur on GPU! - GPUArraysCore.allowscalar(false) - - @test minimum(ϕ) ≈ minimum(ϕ_vals) atol=ε - @test maximum(ϕ) ≈ maximum(ϕ_vals) atol=ε - @test mean(ϕ) ≈ mean(ϕ_vals) atol=2ε - @test minimum(∛, ϕ) ≈ minimum(∛, ϕ_vals) atol=ε - @test maximum(abs, ϕ) ≈ maximum(abs, ϕ_vals) atol=ε - @test mean(abs2, ϕ) ≈ mean(abs2, ϕ) atol=ε - - @test extrema(ϕ) == (minimum(ϕ), maximum(ϕ)) - @test extrema(∛, ϕ) == (minimum(∛, ϕ), maximum(∛, ϕ)) - - for dims in dims_to_test - @test all(isapprox(minimum(ϕ, dims=dims), minimum(ϕ_vals, dims=dims), atol=4ε)) - @test all(isapprox(maximum(ϕ, dims=dims), maximum(ϕ_vals, dims=dims), atol=4ε)) - @test all(isapprox(mean(ϕ, dims=dims), mean(ϕ_vals, dims=dims), atol=4ε)) - - @test all(isapprox(minimum(sin, ϕ, dims=dims), minimum(sin, ϕ_vals, dims=dims), atol=4ε)) - @test all(isapprox(maximum(cos, ϕ, dims=dims), maximum(cos, ϕ_vals, dims=dims), atol=4ε)) - @test all(isapprox(mean(cosh, ϕ, dims=dims), mean(cosh, ϕ_vals, dims=dims), atol=5ε)) - end - end - - return nothing -end - -# Choose a trilinear function so trilinear interpolation can return values that -# are exactly correct. -@inline func(x, y, z) = convert(typeof(x), exp(-1) + 3x - y/7 + z + 2x*y - 3x*z + 4y*z - 5x*y*z) - -function run_field_interpolation_tests(grid) - arch = architecture(grid) - velocities = VelocityFields(grid) - tracers = TracerFields((:c,), grid) - - (u, v, w), c = velocities, tracers.c - - # Maximum expected rounding error is the unit in last place of the maximum value - # of func over the domain of the grid. - - # TODO: remove this allowscalar when `nodes` returns broadcastable object on GPU - xf, yf, zf = nodes(grid, (Face(), Face(), Face()), reshape=true) - f_max = @allowscalar maximum(func.(xf, yf, zf)) - ε_max = eps(f_max) - tolerance = 10 * ε_max - - set!(u, func) - set!(v, func) - set!(w, func) - set!(c, func) - - # Check that interpolating to the field's own grid points returns - # the same value as the field itself. - for f in (u, v, w, c) - loc = Tuple(L() for L in location(f)) - - result = true - @allowscalar begin - for i in size(f, 1), j in size(f, 2), k in size(f, 3) - x, y, z = node(i, j, k, f) - ℑf = interpolate((x, y, z), f, loc, f.grid) - true_value = interior(f, i, j, k)[] - - # If at last one of the points is not approximately equal to the true value, set result to false and break - if !isapprox(ℑf, true_value, atol=tolerance) - result = false - break - end - end - end - @test result - end - - # Check that interpolating between grid points works as expected. - - xs = Array(reshape([0.3, 0.55, 0.73], (3, 1, 1))) - ys = Array(reshape([-π/6, 0, 1+1e-7], (1, 3, 1))) - zs = Array(reshape([-1.3, 1.23, 2.1], (1, 1, 3))) - - X = [(xs[i], ys[j], zs[k]) for i=1:3, j=1:3, k=1:3] - X = on_architecture(arch, X) - - xs = on_architecture(arch, xs) - ys = on_architecture(arch, ys) - zs = on_architecture(arch, zs) - - @allowscalar begin - for f in (u, v, w, c) - loc = Tuple(L() for L in location(f)) - result = true - for i in size(f, 1), j in size(f, 2), k in size(f, 3) - xi, yi, zi = node(i, j, k, f) - ℑf = interpolate((xi, yi, zi), f, loc, f.grid) - true_value = func(xi, yi, zi) - - # If at last one of the points is not approximately equal to the true value, set result to false and break - if !isapprox(ℑf, true_value, atol=tolerance) - result = false - break - end - end - @test result - - # for the next test we first call fill_halo_regions! on the - # original field `f` - # note, that interpolate! will call fill_halo_regions! on - # the interpolated field after the interpolation - fill_halo_regions!(f) - - f_copy = deepcopy(f) - fill!(f_copy, 0) - interpolate!(f_copy, f) - - @test all(interior(f_copy) .≈ interior(f)) - end - end - - @info " Testing the convert functions" - for n in 1:30 - @test convert_to_0_360(- 10.e0^(-n)) > 359 - @test convert_to_0_360(- 10.f0^(-n)) > 359 - @test convert_to_0_360(10.e0^(-n)) < 1 - @test convert_to_0_360(10.f0^(-n)) < 1 - end - - # Generating a random longitude left bound between -1000 and 1000 - λs₀ = rand(1000) .* 2000 .- 1000 - - # Generating a random interpolation longitude - λsᵢ = rand(1000) .* 2000 .- 1000 - - for λ₀ in λs₀, λᵢ in λsᵢ - @test λ₀ ≤ convert_to_λ₀_λ₀_plus360(λᵢ, λ₀) ≤ λ₀ + 360 - end - - # Check interpolation on Windowed fields - wf = ZFaceField(grid; indices=(:, :, grid.Nz+1)) - If = Field{Center, Center, Nothing}(grid) - set!(If, (x, y)-> x * y) - interpolate!(wf, If) - - @allowscalar begin - @test all(interior(wf) .≈ interior(If)) - end - - # interpolation between fields on latitudelongitude grids with different longitudes - grid1 = LatitudeLongitudeGrid(size=(10, 1, 1), longitude=( 0, 360), latitude=(-90, 90), z=(0, 1)) - grid2 = LatitudeLongitudeGrid(size=(10, 1, 1), longitude=( -180, 180), latitude=(-90, 90), z=(0, 1)) - grid3 = LatitudeLongitudeGrid(size=(10, 1, 1), longitude=(-1080, -1080+360), latitude=(-90, 90), z=(0, 1)) - grid4 = LatitudeLongitudeGrid(size=(10, 1, 1), longitude=( 180, 540), latitude=(-90, 90), z=(0, 1)) - - f1 = CenterField(grid1) - f2 = CenterField(grid2) - f3 = CenterField(grid3) - f4 = CenterField(grid4) - - set!(f1, (λ, y, z) -> λ) - fill_halo_regions!(f1) - interpolate!(f2, f1) - interpolate!(f3, f1) - interpolate!(f4, f1) - - @test all(interior(f2) .≈ map(convert_to_0_360, λnodes(grid2, Center()))) - @test all(interior(f3) .≈ map(convert_to_0_360, λnodes(grid3, Center()))) - @test all(interior(f4) .≈ map(convert_to_0_360, λnodes(grid4, Center()))) - - # now interpolate back - fill_halo_regions!(f2) - fill_halo_regions!(f3) - fill_halo_regions!(f4) - - interpolate!(f1, f2) - @test all(interior(f1) .≈ λnodes(grid1, Center())) - - interpolate!(f1, f3) - @test all(interior(f1) .≈ λnodes(grid1, Center())) - - interpolate!(f1, f4) - @test all(interior(f1) .≈ λnodes(grid1, Center())) - - return nothing -end - -function nodes_of_field_views_are_consistent(grid) - # Test with different field types - test_fields = [CenterField(grid), XFaceField(grid), YFaceField(grid), ZFaceField(grid)] - - for field in test_fields - loc = instantiated_location(field) - - # Test various view patterns - test_indices = [ - (2:6, :, :), # x slice - (:, 2:4, :), # y slice - (:, :, 2:3), # z slice - (3:5, 2:4, :), # xy slice - (2:6, :, 2:3), # xz slice - (:, 2:4, 2:3), # yz slice - (3:5, 2:4, 2:3), # xyz slice - ] - - for test_idx in test_indices - # Create field view with these indices - field_view = view(field, test_idx...) - - # Get nodes from the view - view_nodes = nodes(field_view) - - # Get nodes from the original field with the same indices - # This is what should be equivalent to the view_nodes - full_nodes = nodes(field.grid, loc...; indices=test_idx) - - # Test that they are equal - @test view_nodes == full_nodes - - # Also test that the view's indices match what we expect - @test indices(field_view) == test_idx - - # Test that view nodes have sizes consistent with the view indices - for (i, coord_nodes) in enumerate(view_nodes) - if coord_nodes !== nothing && full_nodes[i] !== nothing - @test coord_nodes == full_nodes[i] - end - end - end - end - - return nothing -end - - -##### -##### -##### - -@testset "Fields" begin - @info "Testing Fields..." - - @testset "Field initialization" begin - @info " Testing Field initialization..." - - N = (4, 6, 8) - L = (2π, 3π, 5π) - H = (1, 1, 1) - - for arch in archs, FT in float_types - grid = RectilinearGrid(arch, FT, size=N, extent=L, halo=H, topology=(Periodic, Periodic, Periodic)) - @test correct_field_size(grid, (Center, Center, Center), N[1] + 2 * H[1], N[2] + 2 * H[2], N[3] + 2 * H[3]) - @test correct_field_size(grid, (Face, Center, Center), N[1] + 2 * H[1], N[2] + 2 * H[2], N[3] + 2 * H[3]) - @test correct_field_size(grid, (Center, Face, Center), N[1] + 2 * H[1], N[2] + 2 * H[2], N[3] + 2 * H[3]) - @test correct_field_size(grid, (Center, Center, Face), N[1] + 2 * H[1], N[2] + 2 * H[2], N[3] + 2 * H[3]) - - grid = RectilinearGrid(arch, FT, size=N, extent=L, halo=H, topology=(Periodic, Periodic, Bounded)) - @test correct_field_size(grid, (Center, Center, Center), N[1] + 2 * H[1], N[2] + 2 * H[2], N[3] + 2 * H[3]) - @test correct_field_size(grid, (Face, Center, Center), N[1] + 2 * H[1], N[2] + 2 * H[2], N[3] + 2 * H[3]) - @test correct_field_size(grid, (Center, Face, Center), N[1] + 2 * H[1], N[2] + 2 * H[2], N[3] + 2 * H[3]) - @test correct_field_size(grid, (Center, Center, Face), N[1] + 2 * H[1], N[2] + 2 * H[2], N[3] + 2 * H[3] + 1) - - grid = RectilinearGrid(arch, FT, size=N, extent=L, halo=H, topology=(Periodic, Bounded, Bounded)) - @test correct_field_size(grid, (Center, Center, Center), N[1] + 2 * H[1], N[2] + 2 * H[2], N[3] + 2 * H[3]) - @test correct_field_size(grid, (Face, Center, Center), N[1] + 2 * H[1], N[2] + 2 * H[2], N[3] + 2 * H[3]) - @test correct_field_size(grid, (Center, Face, Center), N[1] + 2 * H[1], N[2] + 1 + 2 * H[2], N[3] + 2 * H[3]) - @test correct_field_size(grid, (Center, Center, Face), N[1] + 2 * H[1], N[2] + 2 * H[2], N[3] + 1 + 2 * H[3]) - - grid = RectilinearGrid(arch, FT, size=N, extent=L, halo=H, topology=(Bounded, Bounded, Bounded)) - @test correct_field_size(grid, (Center, Center, Center), N[1] + 2 * H[1], N[2] + 2 * H[2], N[3] + 2 * H[3]) - @test correct_field_size(grid, (Face, Center, Center), N[1] + 1 + 2 * H[1], N[2] + 2 * H[2], N[3] + 2 * H[3]) - @test correct_field_size(grid, (Center, Face, Center), N[1] + 2 * H[1], N[2] + 1 + 2 * H[2], N[3] + 2 * H[3]) - @test correct_field_size(grid, (Center, Center, Face), N[1] + 2 * H[1], N[2] + 2 * H[2], N[3] + 1 + 2 * H[3]) - - # Reduced fields - @test correct_field_size(grid, (Nothing, Center, Center), 1, N[2] + 2 * H[2], N[3] + 2 * H[3]) - @test correct_field_size(grid, (Nothing, Center, Center), 1, N[2] + 2 * H[2], N[3] + 2 * H[3]) - @test correct_field_size(grid, (Nothing, Face, Center), 1, N[2] + 2 * H[2] + 1, N[3] + 2 * H[3]) - @test correct_field_size(grid, (Nothing, Face, Face), 1, N[2] + 2 * H[2] + 1, N[3] + 2 * H[3] + 1) - @test correct_field_size(grid, (Center, Nothing, Center), N[1] + 2 * H[1], 1, N[3] + 2 * H[3]) - @test correct_field_size(grid, (Center, Nothing, Center), N[1] + 2 * H[1], 1, N[3] + 2 * H[3]) - @test correct_field_size(grid, (Center, Center, Nothing), N[1] + 2 * H[1], N[2] + 2 * H[2], 1) - @test correct_field_size(grid, (Nothing, Nothing, Center), 1, 1, N[3] + 2 * H[3]) - @test correct_field_size(grid, (Center, Nothing, Nothing), N[1] + 2 * H[1], 1, 1) - @test correct_field_size(grid, (Nothing, Nothing, Nothing), 1, 1, 1) - - # "View" fields - for f in [CenterField(grid), XFaceField(grid), YFaceField(grid), ZFaceField(grid)] - - test_indices = [(:, :, :), (1:2, 3:4, 5:6), (1, 1:6, :)] - test_field_sizes = [size(f), (2, 2, 2), (1, 6, size(f, 3))] - test_parent_sizes = [size(parent(f)), (2, 2, 2), (1, 6, size(parent(f), 3))] - - for (t, indices) in enumerate(test_indices) - field_sz = test_field_sizes[t] - parent_sz = test_parent_sizes[t] - f_view = view(f, indices...) - f_sliced = Field(f; indices) - @test size(f_view) == field_sz - @test size(parent(f_view)) == parent_sz - end - end - - grid = RectilinearGrid(arch, FT, size=N, extent=L, halo=H, topology=(Periodic, Periodic, Periodic)) - for side in (:east, :west, :north, :south, :top, :bottom) - for wrong_bc in (ValueBoundaryCondition(0), - FluxBoundaryCondition(0), - GradientBoundaryCondition(0)) - - wrong_kw = Dict(side => wrong_bc) - wrong_bcs = FieldBoundaryConditions(grid, (Center(), Center(), Center()); wrong_kw...) - @test_throws ArgumentError CenterField(grid, boundary_conditions=wrong_bcs) - end - end - - grid = RectilinearGrid(arch, FT, size=N[2:3], extent=L[2:3], halo=H[2:3], topology=(Flat, Periodic, Periodic)) - for side in (:east, :west) - for wrong_bc in (ValueBoundaryCondition(0), - FluxBoundaryCondition(0), - GradientBoundaryCondition(0)) - - wrong_kw = Dict(side => wrong_bc) - wrong_bcs = FieldBoundaryConditions(grid, (Center(), Center(), Center()); wrong_kw...) - @test_throws ArgumentError CenterField(grid, boundary_conditions=wrong_bcs) - end - end - - grid = RectilinearGrid(arch, FT, size=N, extent=L, halo=H, topology=(Periodic, Bounded, Bounded)) - for side in (:east, :west, :north, :south) - for wrong_bc in (ValueBoundaryCondition(0), - FluxBoundaryCondition(0), - GradientBoundaryCondition(0)) - - wrong_kw = Dict(side => wrong_bc) - wrong_bcs = FieldBoundaryConditions(grid, (Center(), Face(), Face()); wrong_kw...) - - @test_throws ArgumentError Field{Center, Face, Face}(grid, boundary_conditions=wrong_bcs) - end - end - - if arch isa GPU - wrong_bcs = FieldBoundaryConditions(grid, (Center(), Center(), Center()), - top=FluxBoundaryCondition(zeros(FT, N[1], N[2]))) - @test_throws ArgumentError CenterField(grid, boundary_conditions=wrong_bcs) - end - end - end - - @testset "Setting fields" begin - @info " Testing field setting..." - - FieldTypes = (CenterField, XFaceField, YFaceField, ZFaceField) - - N = (4, 6, 8) - L = (2π, 3π, 5π) - H = (1, 1, 1) - - int_vals = Any[0, Int8(-1), Int16(2), Int32(-3), Int64(4)] - uint_vals = Any[6, UInt8(7), UInt16(8), UInt32(9), UInt64(10)] - float_vals = Any[0.0, -0.0, 6e-34, 1.0f10] - rational_vals = Any[1//11, -23//7] - other_vals = Any[π] - vals = vcat(int_vals, uint_vals, float_vals, rational_vals, other_vals) - - for arch in archs, FT in float_types - ArrayType = array_type(arch) - grid = RectilinearGrid(arch, FT, size=N, extent=L, topology=(Periodic, Periodic, Bounded)) - - for FieldType in FieldTypes, val in vals - @test correct_field_value_was_set(grid, FieldType, val) - end - - for loc in ((Center, Center, Center), - (Face, Center, Center), - (Center, Face, Center), - (Center, Center, Face), - (Nothing, Center, Center), - (Center, Nothing, Center), - (Center, Center, Nothing), - (Nothing, Nothing, Center), - (Nothing, Nothing, Nothing)) - - field = Field(instantiate(loc), grid) - sz = size(field) - A = rand(FT, sz...) - set!(field, A) - @test @allowscalar field.data[1, 1, 1] == A[1, 1, 1] - end - - Nx = 8 - topo = (Bounded, Bounded, Bounded) - grid = RectilinearGrid(arch, FT, topology=topo, size=(Nx, Nx, Nx), x=(-1, 1), y=(0, 2π), z=(-1, 1)) - - @info " Testing field construction with `field` function..." - - array_data = ones(FT, Nx, Nx, Nx) - f = field((Center, Center, Center), array_data, grid) - @test @allowscalar all(isone, interior(f)) - - # With an OffsetArray or a Field, we point to the same data - offset_data = Oceananigans.Grids.new_data(FT, grid, (Center(), Center(), Center())) - fill!(offset_data, 1) - f = field((Center, Center, Center), offset_data, grid) - @test @allowscalar all(isone, f.data) - @test f.data === offset_data - - field_data = CenterField(grid) - set!(field_data, 1) - f = field((Center, Center, Center), field_data, grid) - @test @allowscalar all(isone, interior(f)) - @test f === field_data - - number_data = FT(1) - f = field((Center, Center, Center), number_data, grid) - @test f.constant == 1 - - function_data = (x, y, z) -> 1 - f = field((Center, Center, Center), function_data, grid) - @test @allowscalar all(isone, interior(f)) - - @info " Testing Field constructors..." - - u = XFaceField(grid) - v = YFaceField(grid) - w = ZFaceField(grid) - c = CenterField(grid) - - f(x, y, z) = exp(x) * sin(y) * tanh(z) - - ϕs = (u, v, w, c) - [set!(ϕ, f) for ϕ in ϕs] - - xu, yu, zu = nodes(u) - xv, yv, zv = nodes(v) - xw, yw, zw = nodes(w) - xc, yc, zc = nodes(c) - - @test @allowscalar u[1, 2, 3] ≈ f(xu[1], yu[2], zu[3]) - @test @allowscalar v[1, 2, 3] ≈ f(xv[1], yv[2], zv[3]) - @test @allowscalar w[1, 2, 3] ≈ f(xw[1], yw[2], zw[3]) - @test @allowscalar c[1, 2, 3] ≈ f(xc[1], yc[2], zc[3]) - - # Test for Field-to-Field setting on same architecture, and cross architecture. - # The behavior depends on halo size: if the halos of two fields are the same, we can - # (easily) copy halo data over. - # Otherwise, we take the easy way out (for now) and only copy interior data. - big_halo = (3, 3, 3) - small_halo = (1, 1, 1) - domain = (; x=(0, 1), y=(0, 1), z=(0, 1)) - sz = (3, 3, 3) - - grid = RectilinearGrid(arch, FT; halo=big_halo, size=sz, domain...) - a = CenterField(grid) - b = CenterField(grid) - parent(a) .= 1 - set!(b, a) - @test parent(b) == parent(a) - - grid_with_smaller_halo = RectilinearGrid(arch, FT; halo=small_halo, size=sz, domain...) - c = CenterField(grid_with_smaller_halo) - set!(c, a) - @test interior(c) == interior(a) - - # Cross-architecture setting should have similar behavior - if arch isa GPU - cpu_grid = RectilinearGrid(CPU(), FT; halo=big_halo, size=sz, domain...) - d = CenterField(cpu_grid) - set!(d, a) - @test parent(d) == Array(parent(a)) - @test d == a - @test a == d - - cpu_grid_with_smaller_halo = RectilinearGrid(CPU(), FT; halo=small_halo, size=sz, domain...) - e = CenterField(cpu_grid_with_smaller_halo) - set!(e, a) - @test e == a - @test a == e - @test Array(interior(e)) == Array(interior((a))) - end - end - end - - @testset "isapprox on Fields" begin - for arch in archs, FT in float_types - # Make sure this doesn't require scalar indexing - GPUArraysCore.allowscalar(false) - - rect_grid = RectilinearGrid(arch, FT; size=(8, 8, 8), x=(0, 1_000), y=(0, 1_000), z=(0, 1_000)) - - H = 100.0 - W = 1000.0 - mountain(x, y) = H * exp(-(x^2 + y^2) / 2W^2) - imm_grid = ImmersedBoundaryGrid(rect_grid, GridFittedBottom(mountain)) - - for grid in (rect_grid, imm_grid) - @info " Testing isapprox on fields [$(typeof(arch)), $FT, $(nameof(typeof(grid)))]..." - u = CenterField(grid) - v = CenterField(grid) - set!(u, 1) - set!(v, 1) - Oceananigans.ImmersedBoundaries.mask_immersed_field!(u, 1) - Oceananigans.ImmersedBoundaries.mask_immersed_field!(v, 2) - # Make sure the two fields are the same - @test isapprox(u, v) - @test isapprox(u, v; rtol=0, atol=0) - - set!(v, FT(1.1)) - @test !isapprox(u, v) - @test isapprox(u, v; rtol=0.1) - @test !isapprox(u, v; atol=2.0) - # norm(u) = √512, norm(v) = 1.1 * √512, difference is 0.1 * √512 ∼ 2.26274, - # we use a slightly larger tolerance to make the check successful. - @test isapprox(u, v; atol=2.26275) - end - end - end - - @testset "Field reductions" begin - for arch in archs, FT in float_types - @info " Testing field reductions [$(typeof(arch)), $FT]..." - N = 8 - topo = (Bounded, Bounded, Bounded) - size = (N, N, N) - y = (0, 2π) - z = (-1, 1) - - x = (-1, 1) - regular_grid = RectilinearGrid(arch, FT; topology=topo, size, x, y, z) - - x = range(-1, stop=1, length=N+1) - variably_spaced_grid = RectilinearGrid(arch, FT; topology=topo, size, x, y, z) - - for (name, grid) in [(:regular_grid => regular_grid), - (:variably_spaced_grid => variably_spaced_grid)] - @info " Testing field reductions on $name..." - run_field_reduction_tests(grid) - end - end - - for arch in archs, FT in float_types - @info " Test reductions on WindowedFields [$(typeof(arch)), $FT]..." - - grid = RectilinearGrid(arch, FT, size=(2, 3, 4), x=(0, 1), y=(0, 1), z=(0, 1)) - c = CenterField(grid) - Random.seed!(42) - set!(c, rand(size(c)...)) - - windowed_c = view(c, :, 2:3, 1:2) - - for fun in (sum, maximum, minimum) - @test fun(c) ≈ fun(interior(c)) - @test fun(windowed_c) ≈ fun(interior(windowed_c)) - end - - @test mean(c) ≈ @allowscalar mean(interior(c)) - @test mean(windowed_c) ≈ @allowscalar mean(interior(windowed_c)) - end - end - - @testset "Unit interpolation" begin - for arch in archs - hu = (-1, 1) - hs = range(-1, 1, length=21) - zu = (-100, 0) - zs = range(-100, 0, length=33) - - for latitude in (hu, hs), longitude in (hu, hs), z in (zu, zs), loc in (Center(), Face()) - @info " Testing interpolation for $(latitude) latitude and longitude, $(z) z on $(typeof(loc))s..." - grid = LatitudeLongitudeGrid(arch; size = (20, 20, 32), longitude, latitude, z, halo = (5, 5, 5)) - - # Test random positions, - # set seed for reproducibility - Random.seed!(1234) - Xs = [(2rand()-1, 2rand()-1, -100rand()) for p in 1:20] - - for X in Xs - (x, y, z) = X - fi = @allowscalar FractionalIndices(X, grid, loc, loc, loc) - - i⁻, i⁺, _ = interpolator(fi.i) - j⁻, j⁺, _ = interpolator(fi.j) - k⁻, k⁺, _ = interpolator(fi.k) - - x⁻ = @allowscalar ξnode(i⁻, j⁻, k⁻, grid, loc, loc, loc) - y⁻ = @allowscalar ηnode(i⁻, j⁻, k⁻, grid, loc, loc, loc) - z⁻ = @allowscalar rnode(i⁻, j⁻, k⁻, grid, loc, loc, loc) - - x⁺ = @allowscalar ξnode(i⁺, j⁺, k⁺, grid, loc, loc, loc) - y⁺ = @allowscalar ηnode(i⁺, j⁺, k⁺, grid, loc, loc, loc) - z⁺ = @allowscalar rnode(i⁺, j⁺, k⁺, grid, loc, loc, loc) - - @test x⁻ ≤ x ≤ x⁺ - @test y⁻ ≤ y ≤ y⁺ - @test z⁻ ≤ z ≤ z⁺ - end - end - end - end - - @testset "Field interpolation" begin - for arch in archs, FT in float_types - @info " Testing field interpolation [$(typeof(arch)), $FT]..." - reg_grid = RectilinearGrid(arch, FT, size=(4, 5, 7), x=(0, 1), y=(-π, π), z=(-5.3, 2.7), halo=(1, 1, 1)) - - # Choose points z points to be rounded values of `reg_grid` z nodes so that interpolation matches tolerance - stretched_grid = RectilinearGrid(arch, - size = (4, 5, 7), - halo = (1, 1, 1), - x = [0.0, 0.26, 0.49, 0.78, 1.0], - y = [-3.1, -1.9, -0.6, 0.6, 1.9, 3.1], - z = [-5.3, -4.2, -3.0, -1.9, -0.7, 0.4, 1.6, 2.7]) - - grids = [reg_grid, stretched_grid] - - for grid in grids - run_field_interpolation_tests(grid) - end - end - end - - @testset "Field utils" begin - @info " Testing field utils..." - - @test has_velocities(()) == false - @test has_velocities((:u,)) == false - @test has_velocities((:u, :v)) == false - @test has_velocities((:u, :v, :w)) == true - - @info " Testing similar(f) for f::Union(Field, ReducedField)..." - - grid = RectilinearGrid(CPU(), size=(1, 1, 1), extent=(1, 1, 1)) - - for X in (Center, Face), Y in (Center, Face), Z in (Center, Face) - for arch in archs - f = Field{X, Y, Z}(grid) - run_similar_field_tests(f) - - for dims in (3, (1, 2), (1, 2, 3)) - loc = reduced_location((X(), Y(), Z()); dims) - f = Field(loc, grid) - run_similar_field_tests(f) - end - end - end - end - - @testset "Views of field views" begin - @info " Testing views of field views..." - - Nx, Ny, Nz = 1, 1, 7 - - FieldTypes = (CenterField, XFaceField, YFaceField, ZFaceField) - ZTopologies = (Periodic, Bounded) - - for arch in archs, FT in float_types, FieldType in FieldTypes, ZTopology in ZTopologies - grid = RectilinearGrid(arch, FT, size=(Nx, Ny, Nz), x=(0, 1), y=(0, 1), z=(0, 1), topology = (Periodic, Periodic, ZTopology)) - Hx, Hy, Hz = halo_size(grid) - - c = FieldType(grid) - set!(c, (x, y, z) -> rand()) - - k_top = total_length(location(c, 3)(), topology(c, 3)(), size(grid, 3)) - - # First test that the regular view is correct - cv = view(c, :, :, 1+1:k_top-1) - @test size(cv) == (Nx, Ny, k_top-2) - @test size(parent(cv)) == (Nx+2Hx, Ny+2Hy, k_top-2) - @allowscalar @test all(cv[i, j, k] == c[i, j, k] for k in 1+1:k_top-1, j in 1:Ny, i in 1:Nx) - - # Now test the views of views - cvv = view(cv, :, :, 1+2:k_top-2) - @test size(cvv) == (Nx, Ny, k_top-4) - @test size(parent(cvv)) == (Nx+2Hx, Ny+2Hy, k_top-4) - @allowscalar @test all(cvv[i, j, k] == cv[i, j, k] for k in 1+2:k_top-2, j in 1:Ny, i in 1:Nx) - - cvvv = view(cvv, :, :, 1+3:k_top-3) - @test size(cvvv) == (1, 1, k_top-6) - @test size(parent(cvvv)) == (Nx+2Hx, Ny+2Hy, k_top-6) - @allowscalar @test all(cvvv[i, j, k] == cvv[i, j, k] for k in 1+3:k_top-3, j in 1:Ny, i in 1:Nx) - - @test_throws ArgumentError view(cv, :, :, 1) - @test_throws ArgumentError view(cv, :, :, k_top) - @test_throws ArgumentError view(cvv, :, :, 1:1+1) - @test_throws ArgumentError view(cvv, :, :, k_top-1:k_top) - @test_throws ArgumentError view(cvvv, :, :, 1:1+2) - @test_throws ArgumentError view(cvvv, :, :, k_top-2:k_top) - - @test_throws BoundsError cv[:, :, 1] - @test_throws BoundsError cv[:, :, k_top] - @test_throws BoundsError cvv[:, :, 1:1+1] - @test_throws BoundsError cvv[:, :, k_top-1:k_top] - @test_throws BoundsError cvvv[:, :, 1:1+2] - @test_throws BoundsError cvvv[:, :, k_top-2:k_top] - end - end - - @testset "Field nodes and view consistency" begin - @info " Testing that nodes() returns indices consistent with view()..." - for arch in archs, FT in float_types - # Test RectilinearGrid - rectilinear_grid = RectilinearGrid(arch, FT, size=(8, 6, 4), extent=(2, 3, 1)) - nodes_of_field_views_are_consistent(rectilinear_grid) - - # Test LatitudeLongitudeGrid - latlon_grid = LatitudeLongitudeGrid(arch, FT, size=(8, 6, 4), longitude = (-180, 180), latitude = (-85, 85), z = (-100, 0)) - nodes_of_field_views_are_consistent(latlon_grid) - - # Test OrthogonalSphericalShellGrid (TripolarGrid) - tripolar_grid = TripolarGrid(arch, FT, size=(8, 6, 4)) - nodes_of_field_views_are_consistent(tripolar_grid) - - # Test Flat topology behavior for RectilinearGrid - flat_rlgrid = RectilinearGrid(arch, FT, size=(), extent=(), topology=(Flat, Flat, Flat)) - c_flat = CenterField(flat_rlgrid) - @test nodes(c_flat) == (nothing, nothing, nothing) - - # Test Flat topology behavior for LatitudeLongitudeGrid - flat_llgrid = LatitudeLongitudeGrid(arch, FT, size=(), topology=(Flat, Flat, Flat)) - c_flat = CenterField(flat_llgrid) - @test nodes(c_flat) == (nothing, nothing, nothing) - end - end -end +include("dependencies_for_runtests.jl") + +using Statistics + +using Oceananigans.Fields: CenterField, ReducedField, has_velocities +using Oceananigans.Fields: VelocityFields, TracerFields, interpolate, interpolate! +using Oceananigans.Fields: reduced_location +using Oceananigans.Fields: FractionalIndices, interpolator, instantiate +using Oceananigans.Fields: convert_to_0_360, convert_to_λ₀_λ₀_plus360 +using Oceananigans.Grids: ξnode, ηnode, rnode +using Oceananigans.Grids: total_length +using Oceananigans.Grids: λnode +using Oceananigans.Grids: RectilinearGrid +using Oceananigans.ImmersedBoundaries: ImmersedBoundaryGrid, GridFittedBottom + +using Random +using GPUArraysCore: @allowscalar + +""" + correct_field_size(grid, FieldType, Tx, Ty, Tz) + +Test that the field initialized by the FieldType constructor on `grid` +has size `(Tx, Ty, Tz)`. +""" +correct_field_size(grid, loc, Tx, Ty, Tz) = size(parent(Field(instantiate(loc), grid))) == (Tx, Ty, Tz) + +function run_similar_field_tests(f) + g = similar(f) + @test typeof(f) == typeof(g) + @test f.grid == g.grid + @test location(f) === location(g) + @test !(f.data === g.data) + return nothing +end + +""" + correct_field_value_was_set(N, L, ftf, val) + +Test that the field initialized by the field type function `ftf` on the grid g +can be correctly filled with the value `val` using the `set!(f::AbstractField, v)` +function. +""" +function correct_field_value_was_set(grid, FieldType, val::Number) + arch = architecture(grid) + f = FieldType(grid) + set!(f, val) + return all(interior(f) .≈ val * on_architecture(arch, ones(size(f)))) +end + +function run_field_reduction_tests(grid) + u = XFaceField(grid) + v = YFaceField(grid) + w = ZFaceField(grid) + c = CenterField(grid) + η = Field{Center, Center, Nothing}(grid) + + f(x, y, z) = 1 + exp(x) * sin(y) * tanh(z) + + ϕs = [u, v, w, c] + [set!(ϕ, f) for ϕ in ϕs] + + z_top = znodes(grid, Face())[end] + set!(η, (x, y) -> f(x, y, z_top)) + push!(ϕs, η) + ϕs = Tuple(ϕs) + + u_vals = f.(nodes(u, reshape=true)...) + v_vals = f.(nodes(v, reshape=true)...) + w_vals = f.(nodes(w, reshape=true)...) + c_vals = f.(nodes(c, reshape=true)...) + η_vals = f.(nodes(η, reshape=true)...) + + # Convert to CuArray if needed. + arch = architecture(grid) + u_vals = on_architecture(arch, u_vals) + v_vals = on_architecture(arch, v_vals) + w_vals = on_architecture(arch, w_vals) + c_vals = on_architecture(arch, c_vals) + η_vals = on_architecture(arch, η_vals) + + ϕs_vals = (u_vals, v_vals, w_vals, c_vals, η_vals) + + dims_to_test = (1, 2, 3, (1, 2), (1, 3), (2, 3), (1, 2, 3)) + + for (ϕ, ϕ_vals) in zip(ϕs, ϕs_vals) + ε = eps(eltype(grid)) * 10 * maximum(maximum.(ϕs_vals)) + @info " Testing field reductions with tolerance $ε..." + + @test @allowscalar all(isapprox.(ϕ, ϕ_vals, atol=ε)) # if this isn't true, reduction tests can't pass + + # Important to make sure no scalar operations occur on GPU! + GPUArraysCore.allowscalar(false) + + @test minimum(ϕ) ≈ minimum(ϕ_vals) atol=ε + @test maximum(ϕ) ≈ maximum(ϕ_vals) atol=ε + @test mean(ϕ) ≈ mean(ϕ_vals) atol=2ε + @test minimum(∛, ϕ) ≈ minimum(∛, ϕ_vals) atol=ε + @test maximum(abs, ϕ) ≈ maximum(abs, ϕ_vals) atol=ε + @test mean(abs2, ϕ) ≈ mean(abs2, ϕ) atol=ε + + @test extrema(ϕ) == (minimum(ϕ), maximum(ϕ)) + @test extrema(∛, ϕ) == (minimum(∛, ϕ), maximum(∛, ϕ)) + + for dims in dims_to_test + @test all(isapprox(minimum(ϕ, dims=dims), minimum(ϕ_vals, dims=dims), atol=4ε)) + @test all(isapprox(maximum(ϕ, dims=dims), maximum(ϕ_vals, dims=dims), atol=4ε)) + @test all(isapprox(mean(ϕ, dims=dims), mean(ϕ_vals, dims=dims), atol=4ε)) + + @test all(isapprox(minimum(sin, ϕ, dims=dims), minimum(sin, ϕ_vals, dims=dims), atol=4ε)) + @test all(isapprox(maximum(cos, ϕ, dims=dims), maximum(cos, ϕ_vals, dims=dims), atol=4ε)) + @test all(isapprox(mean(cosh, ϕ, dims=dims), mean(cosh, ϕ_vals, dims=dims), atol=5ε)) + end + end + + return nothing +end + +# Choose a trilinear function so trilinear interpolation can return values that +# are exactly correct. +@inline func(x, y, z) = convert(typeof(x), exp(-1) + 3x - y/7 + z + 2x*y - 3x*z + 4y*z - 5x*y*z) + +function run_field_interpolation_tests(grid) + arch = architecture(grid) + velocities = VelocityFields(grid) + tracers = TracerFields((:c,), grid) + + (u, v, w), c = velocities, tracers.c + + # Maximum expected rounding error is the unit in last place of the maximum value + # of func over the domain of the grid. + + # TODO: remove this allowscalar when `nodes` returns broadcastable object on GPU + xf, yf, zf = nodes(grid, (Face(), Face(), Face()), reshape=true) + f_max = @allowscalar maximum(func.(xf, yf, zf)) + ε_max = eps(f_max) + tolerance = 10 * ε_max + + set!(u, func) + set!(v, func) + set!(w, func) + set!(c, func) + + # Check that interpolating to the field's own grid points returns + # the same value as the field itself. + for f in (u, v, w, c) + loc = Tuple(L() for L in location(f)) + + result = true + @allowscalar begin + for i in size(f, 1), j in size(f, 2), k in size(f, 3) + x, y, z = node(i, j, k, f) + ℑf = interpolate((x, y, z), f, loc, f.grid) + true_value = interior(f, i, j, k)[] + + # If at last one of the points is not approximately equal to the true value, set result to false and break + if !isapprox(ℑf, true_value, atol=tolerance) + result = false + break + end + end + end + @test result + end + + # Check that interpolating between grid points works as expected. + + xs = Array(reshape([0.3, 0.55, 0.73], (3, 1, 1))) + ys = Array(reshape([-π/6, 0, 1+1e-7], (1, 3, 1))) + zs = Array(reshape([-1.3, 1.23, 2.1], (1, 1, 3))) + + X = [(xs[i], ys[j], zs[k]) for i=1:3, j=1:3, k=1:3] + X = on_architecture(arch, X) + + xs = on_architecture(arch, xs) + ys = on_architecture(arch, ys) + zs = on_architecture(arch, zs) + + @allowscalar begin + for f in (u, v, w, c) + loc = Tuple(L() for L in location(f)) + result = true + for i in size(f, 1), j in size(f, 2), k in size(f, 3) + xi, yi, zi = node(i, j, k, f) + ℑf = interpolate((xi, yi, zi), f, loc, f.grid) + true_value = func(xi, yi, zi) + + # If at last one of the points is not approximately equal to the true value, set result to false and break + if !isapprox(ℑf, true_value, atol=tolerance) + result = false + break + end + end + @test result + + # for the next test we first call fill_halo_regions! on the + # original field `f` + # note, that interpolate! will call fill_halo_regions! on + # the interpolated field after the interpolation + fill_halo_regions!(f) + + f_copy = deepcopy(f) + fill!(f_copy, 0) + interpolate!(f_copy, f) + + @test all(interior(f_copy) .≈ interior(f)) + end + end + + @info " Testing the convert functions" + for n in 1:30 + @test convert_to_0_360(- 10.e0^(-n)) > 359 + @test convert_to_0_360(- 10.f0^(-n)) > 359 + @test convert_to_0_360(10.e0^(-n)) < 1 + @test convert_to_0_360(10.f0^(-n)) < 1 + end + + # Generating a random longitude left bound between -1000 and 1000 + λs₀ = rand(1000) .* 2000 .- 1000 + + # Generating a random interpolation longitude + λsᵢ = rand(1000) .* 2000 .- 1000 + + for λ₀ in λs₀, λᵢ in λsᵢ + @test λ₀ ≤ convert_to_λ₀_λ₀_plus360(λᵢ, λ₀) ≤ λ₀ + 360 + end + + # Check interpolation on Windowed fields + wf = ZFaceField(grid; indices=(:, :, grid.Nz+1)) + If = Field{Center, Center, Nothing}(grid) + set!(If, (x, y)-> x * y) + interpolate!(wf, If) + + @allowscalar begin + @test all(interior(wf) .≈ interior(If)) + end + + # interpolation between fields on latitudelongitude grids with different longitudes + grid1 = LatitudeLongitudeGrid(size=(10, 1, 1), longitude=( 0, 360), latitude=(-90, 90), z=(0, 1)) + grid2 = LatitudeLongitudeGrid(size=(10, 1, 1), longitude=( -180, 180), latitude=(-90, 90), z=(0, 1)) + grid3 = LatitudeLongitudeGrid(size=(10, 1, 1), longitude=(-1080, -1080+360), latitude=(-90, 90), z=(0, 1)) + grid4 = LatitudeLongitudeGrid(size=(10, 1, 1), longitude=( 180, 540), latitude=(-90, 90), z=(0, 1)) + + f1 = CenterField(grid1) + f2 = CenterField(grid2) + f3 = CenterField(grid3) + f4 = CenterField(grid4) + + set!(f1, (λ, y, z) -> λ) + fill_halo_regions!(f1) + interpolate!(f2, f1) + interpolate!(f3, f1) + interpolate!(f4, f1) + + @test all(interior(f2) .≈ map(convert_to_0_360, λnodes(grid2, Center()))) + @test all(interior(f3) .≈ map(convert_to_0_360, λnodes(grid3, Center()))) + @test all(interior(f4) .≈ map(convert_to_0_360, λnodes(grid4, Center()))) + + # now interpolate back + fill_halo_regions!(f2) + fill_halo_regions!(f3) + fill_halo_regions!(f4) + + interpolate!(f1, f2) + @test all(interior(f1) .≈ λnodes(grid1, Center())) + + interpolate!(f1, f3) + @test all(interior(f1) .≈ λnodes(grid1, Center())) + + interpolate!(f1, f4) + @test all(interior(f1) .≈ λnodes(grid1, Center())) + + return nothing +end + +function nodes_of_field_views_are_consistent(grid) + # Test with different field types + test_fields = [CenterField(grid), XFaceField(grid), YFaceField(grid), ZFaceField(grid)] + + for field in test_fields + loc = instantiated_location(field) + + # Test various view patterns + test_indices = [ + (2:6, :, :), # x slice + (:, 2:4, :), # y slice + (:, :, 2:3), # z slice + (3:5, 2:4, :), # xy slice + (2:6, :, 2:3), # xz slice + (:, 2:4, 2:3), # yz slice + (3:5, 2:4, 2:3), # xyz slice + ] + + for test_idx in test_indices + # Create field view with these indices + field_view = view(field, test_idx...) + + # Get nodes from the view + view_nodes = nodes(field_view) + + # Get nodes from the original field with the same indices + # This is what should be equivalent to the view_nodes + full_nodes = nodes(field.grid, loc...; indices=test_idx) + + # Test that they are equal + @test view_nodes == full_nodes + + # Also test that the view's indices match what we expect + @test indices(field_view) == test_idx + + # Test that view nodes have sizes consistent with the view indices + for (i, coord_nodes) in enumerate(view_nodes) + if coord_nodes !== nothing && full_nodes[i] !== nothing + @test coord_nodes == full_nodes[i] + end + end + end + end + + return nothing +end + + +##### +##### +##### + +@testset "Fields" begin + @info "Testing Fields..." + + @testset "Field initialization" begin + @info " Testing Field initialization..." + + N = (4, 6, 8) + L = (2π, 3π, 5π) + H = (1, 1, 1) + + for arch in archs, FT in float_types + grid = RectilinearGrid(arch, FT, size=N, extent=L, halo=H, topology=(Periodic, Periodic, Periodic)) + @test correct_field_size(grid, (Center, Center, Center), N[1] + 2 * H[1], N[2] + 2 * H[2], N[3] + 2 * H[3]) + @test correct_field_size(grid, (Face, Center, Center), N[1] + 2 * H[1], N[2] + 2 * H[2], N[3] + 2 * H[3]) + @test correct_field_size(grid, (Center, Face, Center), N[1] + 2 * H[1], N[2] + 2 * H[2], N[3] + 2 * H[3]) + @test correct_field_size(grid, (Center, Center, Face), N[1] + 2 * H[1], N[2] + 2 * H[2], N[3] + 2 * H[3]) + + grid = RectilinearGrid(arch, FT, size=N, extent=L, halo=H, topology=(Periodic, Periodic, Bounded)) + @test correct_field_size(grid, (Center, Center, Center), N[1] + 2 * H[1], N[2] + 2 * H[2], N[3] + 2 * H[3]) + @test correct_field_size(grid, (Face, Center, Center), N[1] + 2 * H[1], N[2] + 2 * H[2], N[3] + 2 * H[3]) + @test correct_field_size(grid, (Center, Face, Center), N[1] + 2 * H[1], N[2] + 2 * H[2], N[3] + 2 * H[3]) + @test correct_field_size(grid, (Center, Center, Face), N[1] + 2 * H[1], N[2] + 2 * H[2], N[3] + 2 * H[3] + 1) + + grid = RectilinearGrid(arch, FT, size=N, extent=L, halo=H, topology=(Periodic, Bounded, Bounded)) + @test correct_field_size(grid, (Center, Center, Center), N[1] + 2 * H[1], N[2] + 2 * H[2], N[3] + 2 * H[3]) + @test correct_field_size(grid, (Face, Center, Center), N[1] + 2 * H[1], N[2] + 2 * H[2], N[3] + 2 * H[3]) + @test correct_field_size(grid, (Center, Face, Center), N[1] + 2 * H[1], N[2] + 1 + 2 * H[2], N[3] + 2 * H[3]) + @test correct_field_size(grid, (Center, Center, Face), N[1] + 2 * H[1], N[2] + 2 * H[2], N[3] + 1 + 2 * H[3]) + + grid = RectilinearGrid(arch, FT, size=N, extent=L, halo=H, topology=(Bounded, Bounded, Bounded)) + @test correct_field_size(grid, (Center, Center, Center), N[1] + 2 * H[1], N[2] + 2 * H[2], N[3] + 2 * H[3]) + @test correct_field_size(grid, (Face, Center, Center), N[1] + 1 + 2 * H[1], N[2] + 2 * H[2], N[3] + 2 * H[3]) + @test correct_field_size(grid, (Center, Face, Center), N[1] + 2 * H[1], N[2] + 1 + 2 * H[2], N[3] + 2 * H[3]) + @test correct_field_size(grid, (Center, Center, Face), N[1] + 2 * H[1], N[2] + 2 * H[2], N[3] + 1 + 2 * H[3]) + + # Reduced fields + @test correct_field_size(grid, (Nothing, Center, Center), 1, N[2] + 2 * H[2], N[3] + 2 * H[3]) + @test correct_field_size(grid, (Nothing, Center, Center), 1, N[2] + 2 * H[2], N[3] + 2 * H[3]) + @test correct_field_size(grid, (Nothing, Face, Center), 1, N[2] + 2 * H[2] + 1, N[3] + 2 * H[3]) + @test correct_field_size(grid, (Nothing, Face, Face), 1, N[2] + 2 * H[2] + 1, N[3] + 2 * H[3] + 1) + @test correct_field_size(grid, (Center, Nothing, Center), N[1] + 2 * H[1], 1, N[3] + 2 * H[3]) + @test correct_field_size(grid, (Center, Nothing, Center), N[1] + 2 * H[1], 1, N[3] + 2 * H[3]) + @test correct_field_size(grid, (Center, Center, Nothing), N[1] + 2 * H[1], N[2] + 2 * H[2], 1) + @test correct_field_size(grid, (Nothing, Nothing, Center), 1, 1, N[3] + 2 * H[3]) + @test correct_field_size(grid, (Center, Nothing, Nothing), N[1] + 2 * H[1], 1, 1) + @test correct_field_size(grid, (Nothing, Nothing, Nothing), 1, 1, 1) + + # "View" fields + for f in [CenterField(grid), XFaceField(grid), YFaceField(grid), ZFaceField(grid)] + + test_indices = [(:, :, :), (1:2, 3:4, 5:6), (1, 1:6, :)] + test_field_sizes = [size(f), (2, 2, 2), (1, 6, size(f, 3))] + test_parent_sizes = [size(parent(f)), (2, 2, 2), (1, 6, size(parent(f), 3))] + + for (t, indices) in enumerate(test_indices) + field_sz = test_field_sizes[t] + parent_sz = test_parent_sizes[t] + f_view = view(f, indices...) + f_sliced = Field(f; indices) + @test size(f_view) == field_sz + @test size(parent(f_view)) == parent_sz + end + end + + grid = RectilinearGrid(arch, FT, size=N, extent=L, halo=H, topology=(Periodic, Periodic, Periodic)) + for side in (:east, :west, :north, :south, :top, :bottom) + for wrong_bc in (ValueBoundaryCondition(0), + FluxBoundaryCondition(0), + GradientBoundaryCondition(0)) + + wrong_kw = Dict(side => wrong_bc) + wrong_bcs = FieldBoundaryConditions(grid, (Center(), Center(), Center()); wrong_kw...) + @test_throws ArgumentError CenterField(grid, boundary_conditions=wrong_bcs) + end + end + + grid = RectilinearGrid(arch, FT, size=N[2:3], extent=L[2:3], halo=H[2:3], topology=(Flat, Periodic, Periodic)) + for side in (:east, :west) + for wrong_bc in (ValueBoundaryCondition(0), + FluxBoundaryCondition(0), + GradientBoundaryCondition(0)) + + wrong_kw = Dict(side => wrong_bc) + wrong_bcs = FieldBoundaryConditions(grid, (Center(), Center(), Center()); wrong_kw...) + @test_throws ArgumentError CenterField(grid, boundary_conditions=wrong_bcs) + end + end + + grid = RectilinearGrid(arch, FT, size=N, extent=L, halo=H, topology=(Periodic, Bounded, Bounded)) + for side in (:east, :west, :north, :south) + for wrong_bc in (ValueBoundaryCondition(0), + FluxBoundaryCondition(0), + GradientBoundaryCondition(0)) + + wrong_kw = Dict(side => wrong_bc) + wrong_bcs = FieldBoundaryConditions(grid, (Center(), Face(), Face()); wrong_kw...) + + @test_throws ArgumentError Field{Center, Face, Face}(grid, boundary_conditions=wrong_bcs) + end + end + + if arch isa GPU + wrong_bcs = FieldBoundaryConditions(grid, (Center(), Center(), Center()), + top=FluxBoundaryCondition(zeros(FT, N[1], N[2]))) + @test_throws ArgumentError CenterField(grid, boundary_conditions=wrong_bcs) + end + end + end + + @testset "Setting fields" begin + @info " Testing field setting..." + + FieldTypes = (CenterField, XFaceField, YFaceField, ZFaceField) + + N = (4, 6, 8) + L = (2π, 3π, 5π) + H = (1, 1, 1) + + int_vals = Any[0, Int8(-1), Int16(2), Int32(-3), Int64(4)] + uint_vals = Any[6, UInt8(7), UInt16(8), UInt32(9), UInt64(10)] + float_vals = Any[0.0, -0.0, 6e-34, 1.0f10] + rational_vals = Any[1//11, -23//7] + other_vals = Any[π] + vals = vcat(int_vals, uint_vals, float_vals, rational_vals, other_vals) + + for arch in archs, FT in float_types + ArrayType = array_type(arch) + grid = RectilinearGrid(arch, FT, size=N, extent=L, topology=(Periodic, Periodic, Bounded)) + + for FieldType in FieldTypes, val in vals + @test correct_field_value_was_set(grid, FieldType, val) + end + + for loc in ((Center, Center, Center), + (Face, Center, Center), + (Center, Face, Center), + (Center, Center, Face), + (Nothing, Center, Center), + (Center, Nothing, Center), + (Center, Center, Nothing), + (Nothing, Nothing, Center), + (Nothing, Nothing, Nothing)) + + field = Field(instantiate(loc), grid) + sz = size(field) + A = rand(FT, sz...) + set!(field, A) + @test @allowscalar field.data[1, 1, 1] == A[1, 1, 1] + end + + Nx = 8 + topo = (Bounded, Bounded, Bounded) + grid = RectilinearGrid(arch, FT, topology=topo, size=(Nx, Nx, Nx), x=(-1, 1), y=(0, 2π), z=(-1, 1)) + + @info " Testing field construction with `field` function..." + + array_data = ones(FT, Nx, Nx, Nx) + f = field((Center, Center, Center), array_data, grid) + @test @allowscalar all(isone, interior(f)) + + # With an OffsetArray or a Field, we point to the same data + offset_data = Oceananigans.Grids.new_data(FT, grid, (Center(), Center(), Center())) + fill!(offset_data, 1) + f = field((Center, Center, Center), offset_data, grid) + @test @allowscalar all(isone, f.data) + @test f.data === offset_data + + field_data = CenterField(grid) + set!(field_data, 1) + f = field((Center, Center, Center), field_data, grid) + @test @allowscalar all(isone, interior(f)) + @test f === field_data + + number_data = FT(1) + f = field((Center, Center, Center), number_data, grid) + @test f.constant == 1 + + function_data = (x, y, z) -> 1 + f = field((Center, Center, Center), function_data, grid) + @test @allowscalar all(isone, interior(f)) + + @info " Testing Field constructors..." + + u = XFaceField(grid) + v = YFaceField(grid) + w = ZFaceField(grid) + c = CenterField(grid) + + f(x, y, z) = exp(x) * sin(y) * tanh(z) + + ϕs = (u, v, w, c) + [set!(ϕ, f) for ϕ in ϕs] + + xu, yu, zu = nodes(u) + xv, yv, zv = nodes(v) + xw, yw, zw = nodes(w) + xc, yc, zc = nodes(c) + + @test @allowscalar u[1, 2, 3] ≈ f(xu[1], yu[2], zu[3]) + @test @allowscalar v[1, 2, 3] ≈ f(xv[1], yv[2], zv[3]) + @test @allowscalar w[1, 2, 3] ≈ f(xw[1], yw[2], zw[3]) + @test @allowscalar c[1, 2, 3] ≈ f(xc[1], yc[2], zc[3]) + + # Test for Field-to-Field setting on same architecture, and cross architecture. + # The behavior depends on halo size: if the halos of two fields are the same, we can + # (easily) copy halo data over. + # Otherwise, we take the easy way out (for now) and only copy interior data. + big_halo = (3, 3, 3) + small_halo = (1, 1, 1) + domain = (; x=(0, 1), y=(0, 1), z=(0, 1)) + sz = (3, 3, 3) + + grid = RectilinearGrid(arch, FT; halo=big_halo, size=sz, domain...) + a = CenterField(grid) + b = CenterField(grid) + parent(a) .= 1 + set!(b, a) + @test parent(b) == parent(a) + + grid_with_smaller_halo = RectilinearGrid(arch, FT; halo=small_halo, size=sz, domain...) + c = CenterField(grid_with_smaller_halo) + set!(c, a) + @test interior(c) == interior(a) + + # Cross-architecture setting should have similar behavior + if arch isa GPU + cpu_grid = RectilinearGrid(CPU(), FT; halo=big_halo, size=sz, domain...) + d = CenterField(cpu_grid) + set!(d, a) + @test parent(d) == Array(parent(a)) + @test d == a + @test a == d + + cpu_grid_with_smaller_halo = RectilinearGrid(CPU(), FT; halo=small_halo, size=sz, domain...) + e = CenterField(cpu_grid_with_smaller_halo) + set!(e, a) + @test e == a + @test a == e + @test Array(interior(e)) == Array(interior((a))) + end + end + end + + @testset "isapprox on Fields" begin + for arch in archs, FT in float_types + # Make sure this doesn't require scalar indexing + GPUArraysCore.allowscalar(false) + + rect_grid = RectilinearGrid(arch, FT; size=(8, 8, 8), x=(0, 1_000), y=(0, 1_000), z=(0, 1_000)) + + H = 100.0 + W = 1000.0 + mountain(x, y) = H * exp(-(x^2 + y^2) / 2W^2) + imm_grid = ImmersedBoundaryGrid(rect_grid, GridFittedBottom(mountain)) + + for grid in (rect_grid, imm_grid) + @info " Testing isapprox on fields [$(typeof(arch)), $FT, $(nameof(typeof(grid)))]..." + u = CenterField(grid) + v = CenterField(grid) + set!(u, 1) + set!(v, 1) + Oceananigans.ImmersedBoundaries.mask_immersed_field!(u, 1) + Oceananigans.ImmersedBoundaries.mask_immersed_field!(v, 2) + # Make sure the two fields are the same + @test isapprox(u, v) + @test isapprox(u, v; rtol=0, atol=0) + + set!(v, FT(1.1)) + @test !isapprox(u, v) + @test isapprox(u, v; rtol=0.1) + @test !isapprox(u, v; atol=2.0) + # norm(u) = √512, norm(v) = 1.1 * √512, difference is 0.1 * √512 ∼ 2.26274, + # we use a slightly larger tolerance to make the check successful. + @test isapprox(u, v; atol=2.26275) + end + end + end + + @testset "Field reductions" begin + for arch in archs, FT in float_types + @info " Testing field reductions [$(typeof(arch)), $FT]..." + N = 8 + topo = (Bounded, Bounded, Bounded) + size = (N, N, N) + y = (0, 2π) + z = (-1, 1) + + x = (-1, 1) + regular_grid = RectilinearGrid(arch, FT; topology=topo, size, x, y, z) + + x = range(-1, stop=1, length=N+1) + variably_spaced_grid = RectilinearGrid(arch, FT; topology=topo, size, x, y, z) + + for (name, grid) in [(:regular_grid => regular_grid), + (:variably_spaced_grid => variably_spaced_grid)] + @info " Testing field reductions on $name..." + run_field_reduction_tests(grid) + end + end + + for arch in archs, FT in float_types + @info " Test reductions on WindowedFields [$(typeof(arch)), $FT]..." + + grid = RectilinearGrid(arch, FT, size=(2, 3, 4), x=(0, 1), y=(0, 1), z=(0, 1)) + c = CenterField(grid) + Random.seed!(42) + set!(c, rand(size(c)...)) + + windowed_c = view(c, :, 2:3, 1:2) + + for fun in (sum, maximum, minimum) + @test fun(c) ≈ fun(interior(c)) + @test fun(windowed_c) ≈ fun(interior(windowed_c)) + end + + @test mean(c) ≈ @allowscalar mean(interior(c)) + @test mean(windowed_c) ≈ @allowscalar mean(interior(windowed_c)) + end + end + + @testset "Unit interpolation" begin + for arch in archs + hu = (-1, 1) + hs = range(-1, 1, length=21) + zu = (-100, 0) + zs = range(-100, 0, length=33) + + for latitude in (hu, hs), longitude in (hu, hs), z in (zu, zs), loc in (Center(), Face()) + @info " Testing interpolation for $(latitude) latitude and longitude, $(z) z on $(typeof(loc))s..." + grid = LatitudeLongitudeGrid(arch; size = (20, 20, 32), longitude, latitude, z, halo = (5, 5, 5)) + + # Test random positions, + # set seed for reproducibility + Random.seed!(1234) + Xs = [(2rand()-1, 2rand()-1, -100rand()) for p in 1:20] + + for X in Xs + (x, y, z) = X + fi = @allowscalar FractionalIndices(X, grid, loc, loc, loc) + + i⁻, i⁺, _ = interpolator(fi.i) + j⁻, j⁺, _ = interpolator(fi.j) + k⁻, k⁺, _ = interpolator(fi.k) + + x⁻ = @allowscalar ξnode(i⁻, j⁻, k⁻, grid, loc, loc, loc) + y⁻ = @allowscalar ηnode(i⁻, j⁻, k⁻, grid, loc, loc, loc) + z⁻ = @allowscalar rnode(i⁻, j⁻, k⁻, grid, loc, loc, loc) + + x⁺ = @allowscalar ξnode(i⁺, j⁺, k⁺, grid, loc, loc, loc) + y⁺ = @allowscalar ηnode(i⁺, j⁺, k⁺, grid, loc, loc, loc) + z⁺ = @allowscalar rnode(i⁺, j⁺, k⁺, grid, loc, loc, loc) + + @test x⁻ ≤ x ≤ x⁺ + @test y⁻ ≤ y ≤ y⁺ + @test z⁻ ≤ z ≤ z⁺ + end + end + end + end + + @testset "Field interpolation" begin + for arch in archs, FT in float_types + @info " Testing field interpolation [$(typeof(arch)), $FT]..." + reg_grid = RectilinearGrid(arch, FT, size=(4, 5, 7), x=(0, 1), y=(-π, π), z=(-5.3, 2.7), halo=(1, 1, 1)) + + # Choose points z points to be rounded values of `reg_grid` z nodes so that interpolation matches tolerance + stretched_grid = RectilinearGrid(arch, + size = (4, 5, 7), + halo = (1, 1, 1), + x = [0.0, 0.26, 0.49, 0.78, 1.0], + y = [-3.1, -1.9, -0.6, 0.6, 1.9, 3.1], + z = [-5.3, -4.2, -3.0, -1.9, -0.7, 0.4, 1.6, 2.7]) + + grids = [reg_grid, stretched_grid] + + for grid in grids + run_field_interpolation_tests(grid) + end + end + end + + @testset "Field utils" begin + @info " Testing field utils..." + + @test has_velocities(()) == false + @test has_velocities((:u,)) == false + @test has_velocities((:u, :v)) == false + @test has_velocities((:u, :v, :w)) == true + + @info " Testing similar(f) for f::Union(Field, ReducedField)..." + + grid = RectilinearGrid(CPU(), size=(1, 1, 1), extent=(1, 1, 1)) + + for X in (Center, Face), Y in (Center, Face), Z in (Center, Face) + for arch in archs + f = Field{X, Y, Z}(grid) + run_similar_field_tests(f) + + for dims in (3, (1, 2), (1, 2, 3)) + loc = reduced_location((X(), Y(), Z()); dims) + f = Field(loc, grid) + run_similar_field_tests(f) + end + end + end + end + + @testset "Views of field views" begin + @info " Testing views of field views..." + + Nx, Ny, Nz = 1, 1, 7 + + FieldTypes = (CenterField, XFaceField, YFaceField, ZFaceField) + ZTopologies = (Periodic, Bounded) + + for arch in archs, FT in float_types, FieldType in FieldTypes, ZTopology in ZTopologies + grid = RectilinearGrid(arch, FT, size=(Nx, Ny, Nz), x=(0, 1), y=(0, 1), z=(0, 1), topology = (Periodic, Periodic, ZTopology)) + Hx, Hy, Hz = halo_size(grid) + + c = FieldType(grid) + set!(c, (x, y, z) -> rand()) + + k_top = total_length(location(c, 3)(), topology(c, 3)(), size(grid, 3)) + + # First test that the regular view is correct + cv = view(c, :, :, 1+1:k_top-1) + @test size(cv) == (Nx, Ny, k_top-2) + @test size(parent(cv)) == (Nx+2Hx, Ny+2Hy, k_top-2) + @allowscalar @test all(cv[i, j, k] == c[i, j, k] for k in 1+1:k_top-1, j in 1:Ny, i in 1:Nx) + + # Now test the views of views + cvv = view(cv, :, :, 1+2:k_top-2) + @test size(cvv) == (Nx, Ny, k_top-4) + @test size(parent(cvv)) == (Nx+2Hx, Ny+2Hy, k_top-4) + @allowscalar @test all(cvv[i, j, k] == cv[i, j, k] for k in 1+2:k_top-2, j in 1:Ny, i in 1:Nx) + + cvvv = view(cvv, :, :, 1+3:k_top-3) + @test size(cvvv) == (1, 1, k_top-6) + @test size(parent(cvvv)) == (Nx+2Hx, Ny+2Hy, k_top-6) + @allowscalar @test all(cvvv[i, j, k] == cvv[i, j, k] for k in 1+3:k_top-3, j in 1:Ny, i in 1:Nx) + + @test_throws ArgumentError view(cv, :, :, 1) + @test_throws ArgumentError view(cv, :, :, k_top) + @test_throws ArgumentError view(cvv, :, :, 1:1+1) + @test_throws ArgumentError view(cvv, :, :, k_top-1:k_top) + @test_throws ArgumentError view(cvvv, :, :, 1:1+2) + @test_throws ArgumentError view(cvvv, :, :, k_top-2:k_top) + + @test_throws BoundsError cv[:, :, 1] + @test_throws BoundsError cv[:, :, k_top] + @test_throws BoundsError cvv[:, :, 1:1+1] + @test_throws BoundsError cvv[:, :, k_top-1:k_top] + @test_throws BoundsError cvvv[:, :, 1:1+2] + @test_throws BoundsError cvvv[:, :, k_top-2:k_top] + end + end + + @testset "Field nodes and view consistency" begin + @info " Testing that nodes() returns indices consistent with view()..." + for arch in archs, FT in float_types + # Test RectilinearGrid + rectilinear_grid = RectilinearGrid(arch, FT, size=(8, 6, 4), extent=(2, 3, 1)) + nodes_of_field_views_are_consistent(rectilinear_grid) + + # Test LatitudeLongitudeGrid + latlon_grid = LatitudeLongitudeGrid(arch, FT, size=(8, 6, 4), longitude = (-180, 180), latitude = (-85, 85), z = (-100, 0)) + nodes_of_field_views_are_consistent(latlon_grid) + + # Test OrthogonalSphericalShellGrid (TripolarGrid) + tripolar_grid = TripolarGrid(arch, FT, size=(8, 6, 4)) + nodes_of_field_views_are_consistent(tripolar_grid) + + # Test Flat topology behavior for RectilinearGrid + flat_rlgrid = RectilinearGrid(arch, FT, size=(), extent=(), topology=(Flat, Flat, Flat)) + c_flat = CenterField(flat_rlgrid) + @test nodes(c_flat) == (nothing, nothing, nothing) + + # Test Flat topology behavior for LatitudeLongitudeGrid + flat_llgrid = LatitudeLongitudeGrid(arch, FT, size=(), topology=(Flat, Flat, Flat)) + c_flat = CenterField(flat_llgrid) + @test nodes(c_flat) == (nothing, nothing, nothing) + end + end +end diff --git a/test/test_gm_infinite_slope.jl b/test/test_gm_infinite_slope.jl index 8310ee844e5..abd84002b9d 100644 --- a/test/test_gm_infinite_slope.jl +++ b/test/test_gm_infinite_slope.jl @@ -20,7 +20,7 @@ function gm_tracer_remains_finite(arch, FT; skew_flux_formulation, horizontal_di x = (0, 1), z = z_faces, topology = (Bounded, Flat, Bounded)) - + model = HydrostaticFreeSurfaceModel(grid; buoyancy = BuoyancyTracer(), closure = eddy_closure, @@ -35,7 +35,7 @@ function gm_tracer_remains_finite(arch, FT; skew_flux_formulation, horizontal_di y = (0, 1), z = z_faces, topology = (Flat, Bounded, Bounded)) - + model = HydrostaticFreeSurfaceModel(grid; buoyancy = BuoyancyTracer(), closure = eddy_closure, @@ -51,7 +51,7 @@ function gm_tracer_remains_finite(arch, FT; skew_flux_formulation, horizontal_di y = (0, 1), z = z_faces, topology = (Bounded, Bounded, Bounded)) - + model = HydrostaticFreeSurfaceModel(grid; buoyancy = BuoyancyTracer(), closure = eddy_closure, diff --git a/test/test_grids.jl b/test/test_grids.jl index 9e426c21cb4..a7d068d8697 100644 --- a/test/test_grids.jl +++ b/test/test_grids.jl @@ -1043,10 +1043,10 @@ end @testset "Latitude-longitude grid" begin @info " Testing latitude-longitude grid construction errors..." - for FT in float_types - @test_throws ArgumentError LatitudeLongitudeGrid(CPU(), FT, size=(10, 10, 4), longitude=(-180, 180), latitude=(-80, 80), z=[-50.0, -30.0, -20.0, 0.0]) # too few z-faces - @test_throws ArgumentError LatitudeLongitudeGrid(CPU(), FT, size=(10, 10, 4), longitude=(-180, 180), latitude=(-80, 80), z=[-2000.0, -1000.0, -50.0, -30.0, -20.0, 0.0]) # too many z-faces - end + for FT in float_types + @test_throws ArgumentError LatitudeLongitudeGrid(CPU(), FT, size=(10, 10, 4), longitude=(-180, 180), latitude=(-80, 80), z=[-50.0, -30.0, -20.0, 0.0]) # too few z-faces + @test_throws ArgumentError LatitudeLongitudeGrid(CPU(), FT, size=(10, 10, 4), longitude=(-180, 180), latitude=(-80, 80), z=[-2000.0, -1000.0, -50.0, -30.0, -20.0, 0.0]) # too many z-faces + end @info " Testing general latitude-longitude grid..." diff --git a/test/test_implicit_diffusion_diagnostic.jl b/test/test_implicit_diffusion_diagnostic.jl index 5fbecdccbb4..442529fe0b6 100644 --- a/test/test_implicit_diffusion_diagnostic.jl +++ b/test/test_implicit_diffusion_diagnostic.jl @@ -22,8 +22,8 @@ function compute_tracer_dissipation!(sim) Δtc² = sim.model.auxiliary_fields.Δtc² Δtd² = sim.model.auxiliary_fields.Δtd² grid = sim.model.grid - Oceananigans.Utils.launch!(architecture(grid), grid, :xyz, - _compute_dissipation!, + Oceananigans.Utils.launch!(architecture(grid), grid, :xyz, + _compute_dissipation!, Δtc², c⁻, c, Δtd², d⁻, d, grid, sim.Δt) return nothing @@ -60,11 +60,11 @@ function test_implicit_diffusion_diagnostic(arch, dim, timestepper, schedule) Δtc² = CenterField(grid) Δtd² = CenterField(grid) - model = HydrostaticFreeSurfaceModel(grid; - timestepper, - velocities, - tracer_advection, - closure, + model = HydrostaticFreeSurfaceModel(grid; + timestepper, + velocities, + tracer_advection, + closure, tracers=(:c, :d), auxiliary_fields=(; Δtc², c⁻, Δtd², d⁻)) @@ -91,9 +91,9 @@ function test_implicit_diffusion_diagnostic(arch, dim, timestepper, schedule) fd = flatten_dissipation_fields(ϵd) outputs = merge(model.tracers, model.auxiliary_fields, fd, fc) - + # Add both callbacks to the simulation with a schedule - add_callback!(sim, ϵc, schedule) + add_callback!(sim, ϵc, schedule) add_callback!(sim, ϵd, schedule) sim.output_writers[:solution] = JLD2Writer(model, outputs; @@ -106,11 +106,11 @@ function test_implicit_diffusion_diagnostic(arch, dim, timestepper, schedule) run!(sim) - Δtc² = FieldTimeSeries("one_d_simulation_$(dim).jld2", "Δtc²") + Δtc² = FieldTimeSeries("one_d_simulation_$(dim).jld2", "Δtc²") Ac = get_advection_dissipation(Val(dim), :c) Dc = get_diffusion_dissipation(Val(dim), :c) - Δtd² = FieldTimeSeries("one_d_simulation_$(dim).jld2", "Δtd²") + Δtd² = FieldTimeSeries("one_d_simulation_$(dim).jld2", "Δtd²") Ad = get_advection_dissipation(Val(dim), :d) Dd = get_diffusion_dissipation(Val(dim), :d) @@ -118,11 +118,11 @@ function test_implicit_diffusion_diagnostic(arch, dim, timestepper, schedule) ∫closs = [sum(interior(Δtc²[i])) for i in 1:Nt] ∫Ac = [sum(interior(Ac[i])) for i in 1:Nt] - ∫Dc = [sum(interior(Dc[i])) for i in 1:Nt] + ∫Dc = [sum(interior(Dc[i])) for i in 1:Nt] ∫dloss = [sum(interior(Δtd²[i])) for i in 1:Nt] ∫Ad = [sum(interior(Ad[i])) for i in 1:Nt] - ∫Dd = [sum(interior(Dd[i])) for i in 1:Nt] + ∫Dd = [sum(interior(Dd[i])) for i in 1:Nt] for i in 3:Nt-1 @test abs(∫closs[i] - ∫Ac[i] - ∫Dc[i]) < 2e-13 # Arbitrary tolerance, not exactly machine precision diff --git a/test/test_init.jl b/test/test_init.jl index 48a656d980e..8950d1f6e45 100644 --- a/test/test_init.jl +++ b/test/test_init.jl @@ -17,4 +17,3 @@ try CUDA.precompile_runtime() @root CUDA.versioninfo() catch; end - diff --git a/test/test_internal_wave_dynamics.jl b/test/test_internal_wave_dynamics.jl index 9de825bfd98..02d772aa813 100644 --- a/test/test_internal_wave_dynamics.jl +++ b/test/test_internal_wave_dynamics.jl @@ -85,4 +85,3 @@ function internal_wave_dynamics_test(model, solution, Δt) return nothing end - diff --git a/test/test_makie_ext.jl b/test/test_makie_ext.jl index ae15095ff12..7f4ed982c4d 100644 --- a/test/test_makie_ext.jl +++ b/test/test_makie_ext.jl @@ -142,7 +142,7 @@ sequential_data(sz::NTuple{N, Int}) where N = reshape(Float64.(1:prod(sz)), sz.. fig = CairoMakie.Figure() ax = CairoMakie.Axis3(fig[1, 1]; aspect=:data) - + # This should work via the extension plt = surface!(ax, T; colormap=:viridis) @test plt isa CairoMakie.Surface @@ -155,7 +155,7 @@ sequential_data(sz::NTuple{N, Int}) where N = reshape(Float64.(1:prod(sz)), sz.. fig = CairoMakie.Figure() ax = CairoMakie.Axis3(fig[1, 1]; aspect=:data) - + plt = surface!(ax, T; colormap=:viridis) @test plt isa CairoMakie.Surface end @@ -176,7 +176,7 @@ sequential_data(sz::NTuple{N, Int}) where N = reshape(Float64.(1:prod(sz)), sz.. fig = CairoMakie.Figure() ax = CairoMakie.Axis3(fig[1, 1]; aspect=:data) - + # This should work with Observable fields plt = surface!(ax, T_obs; colormap=:viridis) @test plt isa CairoMakie.Surface @@ -199,7 +199,7 @@ sequential_data(sz::NTuple{N, Int}) where N = reshape(Float64.(1:prod(sz)), sz.. fig = CairoMakie.Figure() ax = CairoMakie.Axis3(fig[1, 1]; aspect=:data) - + plt = geo_surface!(ax, T; colormap=:thermal) @test plt isa CairoMakie.Surface end diff --git a/test/test_metal.jl b/test/test_metal.jl index 92ce477706a..1bbeb7dae5a 100644 --- a/test/test_metal.jl +++ b/test/test_metal.jl @@ -39,4 +39,3 @@ Oceananigans.defaults.FloatType = Float32 @test iteration(simulation) == 3 @test time(simulation) == 3minutes end - diff --git a/test/test_mpi_tripolar.jl b/test/test_mpi_tripolar.jl index 73403607e7d..35704b39366 100644 --- a/test/test_mpi_tripolar.jl +++ b/test/test_mpi_tripolar.jl @@ -258,4 +258,4 @@ run_large_pencil_distributed_grid = """ @test all(vs .≈ vp) @test all(cs .≈ cp) @test all(ηs .≈ ηp) -end \ No newline at end of file +end diff --git a/test/test_multi_region_unit.jl b/test/test_multi_region_unit.jl index 5c409e295a4..478379479fe 100644 --- a/test/test_multi_region_unit.jl +++ b/test/test_multi_region_unit.jl @@ -76,4 +76,3 @@ end end end end - diff --git a/test/test_operators.jl b/test/test_operators.jl index ae3dc470706..80f17b0acac 100644 --- a/test/test_operators.jl +++ b/test/test_operators.jl @@ -93,11 +93,11 @@ function test_function_differentiation(T=Float64) grid = RectilinearGrid(CPU(), T; size=(3, 3, 3), x=stretched_f, y=stretched_f, z=stretched_f, topology = (Bounded, Bounded, Bounded)) ∂x_f(i, j, k) = 1 / df(i) * (ϕ²[i, j, k] - ϕ²[i-1, j, k]) - ∂x_c(i, j, k) = 1 / dc(i) * (ϕ²[i+1, j, k] - ϕ²[i, j, k]) + ∂x_c(i, j, k) = 1 / dc(i) * (ϕ²[i+1, j, k] - ϕ²[i, j, k]) ∂y_f(i, j, k) = 1 / df(j) * (ϕ²[i, j, k] - ϕ²[i, j-1, k]) - ∂y_c(i, j, k) = 1 / dc(j) * (ϕ²[i, j+1, k] - ϕ²[i, j, k]) + ∂y_c(i, j, k) = 1 / dc(j) * (ϕ²[i, j+1, k] - ϕ²[i, j, k]) ∂z_f(i, j, k) = 1 / df(k) * (ϕ²[i, j, k] - ϕ²[i, j, k-1]) - ∂z_c(i, j, k) = 1 / dc(k) * (ϕ²[i, j, k+1] - ϕ²[i, j, k]) + ∂z_c(i, j, k) = 1 / dc(k) * (ϕ²[i, j, k+1] - ϕ²[i, j, k]) for ∂x in (∂xᶜᶜᶜ, ∂xᶜᶜᶠ, ∂xᶜᶠᶜ, ∂xᶜᶠᶠ) @test ∂x(2, 2, 2, grid, f, ϕ) == ∂x_c(2, 2, 2) diff --git a/test/test_orthogonal_spherical_shell_time_stepping.jl b/test/test_orthogonal_spherical_shell_time_stepping.jl index 9aa8b760166..d020cf9a1f6 100644 --- a/test/test_orthogonal_spherical_shell_time_stepping.jl +++ b/test/test_orthogonal_spherical_shell_time_stepping.jl @@ -50,4 +50,3 @@ using Oceananigans.OrthogonalSphericalShellGrids: RotatedLatitudeLongitudeGrid @test interior(m1.velocities.u) == interior(m2.velocities.u) @test interior(m1.velocities.v) == interior(m2.velocities.v) end - diff --git a/test/test_reactant.jl b/test/test_reactant.jl index e5dd3e201fb..3cc16222c56 100644 --- a/test/test_reactant.jl +++ b/test/test_reactant.jl @@ -14,50 +14,50 @@ function simple_tendency!(model) grid = model.grid arch = grid.architecture Oceananigans.Utils.launch!( - arch, - grid, - :xyz, - _simple_tendency_kernel!, - model.timestepper.Gⁿ.u, - grid, - model.advection.momentum, - model.velocities) + arch, + grid, + :xyz, + _simple_tendency_kernel!, + model.timestepper.Gⁿ.u, + grid, + model.advection.momentum, + model.velocities) return nothing end @testset "Gu kernel" begin - Nx, Ny, Nz = (10, 10, 10) # number of cells - halo = (7, 7, 7) - longitude = (0, 4) - stretched_longitude = [0, 0.1, 0.2, 0.3, 0.4, 0.6, 1.3, 2.5, 2.6, 3.5, 4.0] - latitude = (0, 4) - z = (-1, 0) - lat_lon_kw = (; size=(Nx, Ny, Nz), halo, longitude, latitude, z) - hydrostatic_model_kw = (; momentum_advection=VectorInvariant(), free_surface=ExplicitFreeSurface()) - - arch = Oceananigans.Architectures.ReactantState() - grid = LatitudeLongitudeGrid(arch; lat_lon_kw...) - model = HydrostaticFreeSurfaceModel(grid; hydrostatic_model_kw...) - - ui = randn(size(model.velocities.u)...) - vi = randn(size(model.velocities.v)...) - set!(model, u=ui, v=vi) - - @jit simple_tendency!(model) - - Gu = model.timestepper.Gⁿ.u - Gv = model.timestepper.Gⁿ.v - Gui = Array(interior(Gu)) - Gvi = Array(interior(Gv)) - - carch = Oceananigans.Architectures.ReactantState() - cgrid = LatitudeLongitudeGrid(carch; lat_lon_kw...) - cmodel = HydrostaticFreeSurfaceModel(cgrid; hydrostatic_model_kw...) - - set!(cmodel, u=ui, v=vi) - - simple_tendency!(cmodel) - @test all(Array(interior(model.timestepper.Gⁿ.u)) .≈ Array(interior(cmodel.timestepper.Gⁿ.u))) + Nx, Ny, Nz = (10, 10, 10) # number of cells + halo = (7, 7, 7) + longitude = (0, 4) + stretched_longitude = [0, 0.1, 0.2, 0.3, 0.4, 0.6, 1.3, 2.5, 2.6, 3.5, 4.0] + latitude = (0, 4) + z = (-1, 0) + lat_lon_kw = (; size=(Nx, Ny, Nz), halo, longitude, latitude, z) + hydrostatic_model_kw = (; momentum_advection=VectorInvariant(), free_surface=ExplicitFreeSurface()) + + arch = Oceananigans.Architectures.ReactantState() + grid = LatitudeLongitudeGrid(arch; lat_lon_kw...) + model = HydrostaticFreeSurfaceModel(grid; hydrostatic_model_kw...) + + ui = randn(size(model.velocities.u)...) + vi = randn(size(model.velocities.v)...) + set!(model, u=ui, v=vi) + + @jit simple_tendency!(model) + + Gu = model.timestepper.Gⁿ.u + Gv = model.timestepper.Gⁿ.v + Gui = Array(interior(Gu)) + Gvi = Array(interior(Gv)) + + carch = Oceananigans.Architectures.ReactantState() + cgrid = LatitudeLongitudeGrid(carch; lat_lon_kw...) + cmodel = HydrostaticFreeSurfaceModel(cgrid; hydrostatic_model_kw...) + + set!(cmodel, u=ui, v=vi) + + simple_tendency!(cmodel) + @test all(Array(interior(model.timestepper.Gⁿ.u)) .≈ Array(interior(cmodel.timestepper.Gⁿ.u))) end ridge(λ, φ) = 0.1 * exp((λ - 2)^2 / 2) @@ -264,4 +264,4 @@ using Oceananigans.OutputReaders: cpu_interpolating_time_indices # Test I can index into a Reactant FieldTimeSeries @test fts[5] isa Field -end \ No newline at end of file +end diff --git a/test/test_reactant_latitude_longitude_grid.jl b/test/test_reactant_latitude_longitude_grid.jl index 272f7339278..042ed2ccdf0 100644 --- a/test/test_reactant_latitude_longitude_grid.jl +++ b/test/test_reactant_latitude_longitude_grid.jl @@ -77,4 +77,3 @@ include("reactant_test_utils.jl") test_reactant_model_correctness(LatitudeLongitudeGrid, HydrostaticFreeSurfaceModel, lat_lon_kw, hydrostatic_model_kw) =# end - diff --git a/test/test_schedules.jl b/test/test_schedules.jl index 54a0160222c..d6226b73ad6 100644 --- a/test/test_schedules.jl +++ b/test/test_schedules.jl @@ -93,4 +93,3 @@ using Oceananigans: initialize! st = SpecifiedTimes(2.5) @test 0.4 ≈ schedule_aligned_time_step(st, fake_clock, Inf) end - diff --git a/test/test_sharded_lat_lon.jl b/test/test_sharded_lat_lon.jl index b85ad0d5bcd..a93109d223f 100644 --- a/test/test_sharded_lat_lon.jl +++ b/test/test_sharded_lat_lon.jl @@ -113,4 +113,3 @@ run_pencil_distributed_grid = """ @test all(cs .≈ cp) @test all(ηs .≈ ηp) end - diff --git a/test/test_sharded_tripolar.jl b/test/test_sharded_tripolar.jl index f6591270cd2..7e883d5d07e 100644 --- a/test/test_sharded_tripolar.jl +++ b/test/test_sharded_tripolar.jl @@ -74,4 +74,4 @@ run_pencil_distributed_grid = """ @test all(vs .≈ vp) @test all(cs .≈ cp) @test all(ηs .≈ ηp) -end \ No newline at end of file +end diff --git a/test/test_stokes_drift.jl b/test/test_stokes_drift.jl index df66cc1e36e..8b4769acd10 100644 --- a/test/test_stokes_drift.jl +++ b/test/test_stokes_drift.jl @@ -57,4 +57,3 @@ end end end end - diff --git a/test/test_zstar_coordinate.jl b/test/test_zstar_coordinate.jl index d304366b3b2..d1eed7d65a5 100644 --- a/test/test_zstar_coordinate.jl +++ b/test/test_zstar_coordinate.jl @@ -220,7 +220,7 @@ end model = HydrostaticFreeSurfaceModel(grid; free_surface, tracers = (:b, :c, :constant), - buoyancy = BuoyancyTracer(), + buoyancy = BuoyancyTracer(), vertical_coordinate = ZStarCoordinate(grid)) bᵢ(x, y, z) = x < grid.Lx / 2 ? 0.06 : 0.01 diff --git a/validation/advection/gaussian_bump_advection.jl b/validation/advection/gaussian_bump_advection.jl index f8f315c9adb..1a9215e6225 100644 --- a/validation/advection/gaussian_bump_advection.jl +++ b/validation/advection/gaussian_bump_advection.jl @@ -110,4 +110,4 @@ progress(sim) = @info "time $(prettytime(sim.model.clock.time)), maximum u: $(ma simulation.callbacks[:progress] = Callback(progress, IterationInterval(20)) -run!(simulation) \ No newline at end of file +run!(simulation) diff --git a/validation/advection/plot_rates_convergence_advection.jl b/validation/advection/plot_rates_convergence_advection.jl index f9e6898f50f..0e14a51155e 100644 --- a/validation/advection/plot_rates_convergence_advection.jl +++ b/validation/advection/plot_rates_convergence_advection.jl @@ -168,4 +168,3 @@ plt = plot_solutions!( pnorm, ROC) savefig(plt, "convergence_rates") - diff --git a/validation/bickley_jet/bickley_utils.jl b/validation/bickley_jet/bickley_utils.jl index ffed5272163..08864f2558b 100644 --- a/validation/bickley_jet/bickley_utils.jl +++ b/validation/bickley_jet/bickley_utils.jl @@ -40,4 +40,3 @@ end bickley_grid(; Nh, arch=CPU(), halo=(3, 3, 3)) = RectilinearGrid(arch; size=(Nh, Nh, 1), halo, x = (-2π, 2π), y=(-2π, 2π), z=(0, 1), topology = (Periodic, Periodic, Bounded)) - diff --git a/validation/biogeochemistry/sediment_entrainment.jl b/validation/biogeochemistry/sediment_entrainment.jl index 0fb43e7d8f1..5ad5c36e0df 100644 --- a/validation/biogeochemistry/sediment_entrainment.jl +++ b/validation/biogeochemistry/sediment_entrainment.jl @@ -149,4 +149,3 @@ display(fig) record(fig, "sediment_entrainment.mp4", 1:Nt, framerate=160) do nn n[] = nn end - diff --git a/validation/biogeochemistry/settling_diffusion.jl b/validation/biogeochemistry/settling_diffusion.jl index 49658130c18..a4f5973b042 100644 --- a/validation/biogeochemistry/settling_diffusion.jl +++ b/validation/biogeochemistry/settling_diffusion.jl @@ -71,4 +71,3 @@ record(fig, "settling_diffusion.mp4", 1:100, framerate=24) do nn run!(simulation) update_plot!(simulation) end - diff --git a/validation/biogeochemistry/simple_plankton_continuous_form_biogeochemistry.jl b/validation/biogeochemistry/simple_plankton_continuous_form_biogeochemistry.jl index 0ba621ab8a4..1caaf434d46 100644 --- a/validation/biogeochemistry/simple_plankton_continuous_form_biogeochemistry.jl +++ b/validation/biogeochemistry/simple_plankton_continuous_form_biogeochemistry.jl @@ -178,4 +178,3 @@ simulation.output_writers[:simple_output] = JLD2Writer(model, outputs; filename, overwrite_existing = true) run!(simulation) - diff --git a/validation/biogeochemistry/sinking_tracer.jl b/validation/biogeochemistry/sinking_tracer.jl index 695aef9fd9c..5b25a44b64f 100644 --- a/validation/biogeochemistry/sinking_tracer.jl +++ b/validation/biogeochemistry/sinking_tracer.jl @@ -40,4 +40,3 @@ simulation.callbacks[:plot] = Callback(update!, IterationInterval(100)) display(fig) run!(simulation) - diff --git a/validation/biogeochemistry/two_reacting_tracers.jl b/validation/biogeochemistry/two_reacting_tracers.jl index 79d2f8b2cd1..13ed54524a9 100644 --- a/validation/biogeochemistry/two_reacting_tracers.jl +++ b/validation/biogeochemistry/two_reacting_tracers.jl @@ -56,4 +56,3 @@ record(fig, "tracer_reactions.mp4", 1:100, framerate=24) do _ run!(simulation) update_plot!(simulation) end - diff --git a/validation/convergence_tests/analyze_forced_fixed_slip.jl b/validation/convergence_tests/analyze_forced_fixed_slip.jl index c8db525081d..3b89638e255 100644 --- a/validation/convergence_tests/analyze_forced_fixed_slip.jl +++ b/validation/convergence_tests/analyze_forced_fixed_slip.jl @@ -91,4 +91,3 @@ for (label, error) in zip(labels, errorses) test_rate_of_convergence(L₁[p], Nx[p], expected=-2.0, atol=0.05, name=name * " L₁") test_rate_of_convergence(L∞[p], Nx[p], expected=-2.0, atol=0.10, name=name * " L∞") end - diff --git a/validation/coriolis/enceladus_convection.jl b/validation/coriolis/enceladus_convection.jl index a2ac07cc081..8c782789048 100644 --- a/validation/coriolis/enceladus_convection.jl +++ b/validation/coriolis/enceladus_convection.jl @@ -104,7 +104,7 @@ for T in [0.0] α_Z = thermal_expansion.(T, S, Z_r, Ref(eos)) β_Z = haline_contraction.(T, S, Z_r, Ref(eos)) lines!(ax3, Z_r ./ 1e3, α_Z ./ β_Z, linewidth = 2, label = "T = $T °C, S = $S psu") - + ρ_Z = ρ.(T, S, Z_r, Ref(eos)) lines!(ax4, Z_r ./ 1e3, ρ_Z, linewidth = 2, label = "T = $T °C, S = $S psu") end diff --git a/validation/distributed_simulations/distributed_scaling/run_tests.sh b/validation/distributed_simulations/distributed_scaling/run_tests.sh index 9c60ba75e86..fd9239e57de 100755 --- a/validation/distributed_simulations/distributed_scaling/run_tests.sh +++ b/validation/distributed_simulations/distributed_scaling/run_tests.sh @@ -14,13 +14,13 @@ # 3) If the system is equipped with nsys profiler it is possible to enable a trace with PROFILE_TRACE=1 (in this file) # # 4) Oceananigans is instantiated with the correct MPI build: -# (run these lines in a gpu node substituting modules and paths) -# $ module load my_cuda_module -# $ module load my_cuda_aware_mpi_module -# $ export JULIA_DEPOT_PATH="/path/to/depot" -# $ export JULIA="/path/to/julia" +# (run these lines in a gpu node substituting modules and paths) +# $ module load my_cuda_module +# $ module load my_cuda_aware_mpi_module +# $ export JULIA_DEPOT_PATH="/path/to/depot" +# $ export JULIA="/path/to/julia" # $ $JULIA --check-bounds=no -e 'using Pkg; Pkg.add("MPIPreferences");' -# $ $JULIA --project --check-bounds=no -e 'using MPIPreferences; MPIPreferences.use_system_binaries()' +# $ $JULIA --project --check-bounds=no -e 'using MPIPreferences; MPIPreferences.use_system_binaries()' # $ $JULIA --project --check-bounds=no -e 'using Pkg; Pkg.build("MPI")' # $ $JULIA --project --check-bounds=no -e 'using Pkg; Pkg.instantiate()' # @@ -54,52 +54,52 @@ export SCALING=weak for RX in 1 2 4 8 16 32 64; do for RY in 1 2 4 8 16 32 64; do - export RX + export RX export RY - if test $SIMULATION = "hydrostatic"; then - if test $SCALING = "weak"; then - # Grid size for Weak scaling tests (Hydrostatic) - export NX=$((1440 * RX)) - export NY=$((600 * RY)) - export NZ=100 - else - # Grid size for Strong scaling tests (Hydrostatic) - export NX=1440 - export NY=600 - export NZ=100 - fi - else - if test $SCALING = "weak"; then - # Grid size for Weak scaling tests (Nonhydrostatic) - export NX=$((512 * RX)) - export NY=$((512 * RY)) - export NZ=256 - else - # Grid size for Strong scaling tests (Nonhydrostatic) - export NX=512 - export NY=512 - export NZ=256 - fi - fi + if test $SIMULATION = "hydrostatic"; then + if test $SCALING = "weak"; then + # Grid size for Weak scaling tests (Hydrostatic) + export NX=$((1440 * RX)) + export NY=$((600 * RY)) + export NZ=100 + else + # Grid size for Strong scaling tests (Hydrostatic) + export NX=1440 + export NY=600 + export NZ=100 + fi + else + if test $SCALING = "weak"; then + # Grid size for Weak scaling tests (Nonhydrostatic) + export NX=$((512 * RX)) + export NY=$((512 * RY)) + export NZ=256 + else + # Grid size for Strong scaling tests (Nonhydrostatic) + export NX=512 + export NY=512 + export NZ=256 + fi + fi - RANKS=$((RX * RY)) + RANKS=$((RX * RY)) - export NNODES=$((RANKS / NGPUS_PER_NODE)) - export NTASKS=$NGPUS_PER_NODE + export NNODES=$((RANKS / NGPUS_PER_NODE)) + export NTASKS=$NGPUS_PER_NODE - echo "" - echo "(RX, RY) = $RX, $RY" - echo "(NX, NY) = $NX, $NY" - echo "(NNODES, NTASKS) = $NNODES, $NTASKS" + echo "" + echo "(RX, RY) = $RX, $RY" + echo "(NX, NY) = $NX, $NY" + echo "(NNODES, NTASKS) = $NNODES, $NTASKS" - # ====================================================== # - # ================== RUN SCALING TEST ================== # - # ====================================================== # + # ====================================================== # + # ================== RUN SCALING TEST ================== # + # ====================================================== # - sbatch -N ${NNODES} --gres=gpu:${NTASKS} --ntasks-per-node=${NTASKS} job_script.sh + sbatch -N ${NNODES} --gres=gpu:${NTASKS} --ntasks-per-node=${NTASKS} job_script.sh - # Use qsub on PBS systems!!! - # qsub pbs_job_script.sh + # Use qsub on PBS systems!!! + # qsub pbs_job_script.sh done done diff --git a/validation/field_time_series_boundary_conditions/generate_input_data.jl b/validation/field_time_series_boundary_conditions/generate_input_data.jl index 241191a349d..d51b52b5847 100644 --- a/validation/field_time_series_boundary_conditions/generate_input_data.jl +++ b/validation/field_time_series_boundary_conditions/generate_input_data.jl @@ -36,4 +36,4 @@ function generate_input_data!(grid, times, filename; ρ_ocean = 1000, cp_ocean = set!(Qˢ, Qˢ_tmp, t) set!(τˣ, τˣ_tmp, t) end -end \ No newline at end of file +end diff --git a/validation/float_precision_tests/hydrostatic_two_dimensional_turbulence.jl b/validation/float_precision_tests/hydrostatic_two_dimensional_turbulence.jl index a3c536b61d4..58dc48afc7d 100644 --- a/validation/float_precision_tests/hydrostatic_two_dimensional_turbulence.jl +++ b/validation/float_precision_tests/hydrostatic_two_dimensional_turbulence.jl @@ -113,4 +113,3 @@ catch err end MPI.Barrier(arch.communicator) - diff --git a/validation/float_precision_tests/two_dimensional_turbulence.jl b/validation/float_precision_tests/two_dimensional_turbulence.jl index fbe80d0f9d3..18b46599db5 100644 --- a/validation/float_precision_tests/two_dimensional_turbulence.jl +++ b/validation/float_precision_tests/two_dimensional_turbulence.jl @@ -50,4 +50,3 @@ ow = JLD2Writer(model, outputs; simulation.output_writers[:jld2] = ow run!(simulation) - diff --git a/validation/immersed_boundaries/cylinder_flow_with_tracer.jl b/validation/immersed_boundaries/cylinder_flow_with_tracer.jl index de109f80d3b..3523ca16043 100644 --- a/validation/immersed_boundaries/cylinder_flow_with_tracer.jl +++ b/validation/immersed_boundaries/cylinder_flow_with_tracer.jl @@ -379,5 +379,3 @@ advection = Centered() experiment_name = run_cylinder_steadystate(Nh = 350, advection = advection, radius = R, stop_time = 100, ν = nu) visualize_cylinder_steadystate(experiment_name) analyze_cylinder_steadystate(experiment_name) - - diff --git a/validation/immersed_boundaries/flow_over_hills.jl b/validation/immersed_boundaries/flow_over_hills.jl index 195cdb48144..14dd7578f90 100644 --- a/validation/immersed_boundaries/flow_over_hills.jl +++ b/validation/immersed_boundaries/flow_over_hills.jl @@ -216,4 +216,3 @@ moviename = @sprintf("flow_over_hills_%dd_h%d.mp4", Nx, 10h) record(fig, moviename, 1:Nt, framerate=24) do nn n[] = nn end - diff --git a/validation/immersed_boundaries/mask_binary_operations.jl b/validation/immersed_boundaries/mask_binary_operations.jl index fda92fb7d01..5dc8f765ce5 100644 --- a/validation/immersed_boundaries/mask_binary_operations.jl +++ b/validation/immersed_boundaries/mask_binary_operations.jl @@ -40,4 +40,3 @@ for (n, f) in enumerate([Σc, Πc, ratio_c, pow_c]) end fig - diff --git a/validation/immersed_boundaries/partial_cell_grid_analysis.jl b/validation/immersed_boundaries/partial_cell_grid_analysis.jl index 577dfd34cd4..ef2d74b768b 100644 --- a/validation/immersed_boundaries/partial_cell_grid_analysis.jl +++ b/validation/immersed_boundaries/partial_cell_grid_analysis.jl @@ -46,4 +46,3 @@ Colorbar(fig[1, 2], hmc) hmf = heatmap!(ax_f, fractional_Δzᶠ) Colorbar(fig[2, 2], hmf) display(fig) - diff --git a/validation/immersed_boundaries/resting_stratified_bumpy_ocean.jl b/validation/immersed_boundaries/resting_stratified_bumpy_ocean.jl index 40fd8469544..df0693f8bc4 100644 --- a/validation/immersed_boundaries/resting_stratified_bumpy_ocean.jl +++ b/validation/immersed_boundaries/resting_stratified_bumpy_ocean.jl @@ -112,4 +112,3 @@ hmvd = heatmap!(ax_vd, Δv) =# display(fig) - diff --git a/validation/immersed_boundaries/staircase_convection.jl b/validation/immersed_boundaries/staircase_convection.jl index 84572c57a8c..0bee1b26ea0 100644 --- a/validation/immersed_boundaries/staircase_convection.jl +++ b/validation/immersed_boundaries/staircase_convection.jl @@ -30,7 +30,7 @@ function setup_grid(Nx, Ny, Nz, arch) zs[2:end-1] .+= randn(length(zs[2:end-1])) * (1 / Nz) / 10 grid = RectilinearGrid(arch, Float64, - size = (Nx, Ny, Nz), + size = (Nx, Ny, Nz), halo = (4, 4, 4), x = (0, 1), y = (0, 1), @@ -38,7 +38,7 @@ function setup_grid(Nx, Ny, Nz, arch) # z = zs, topology = (Bounded, Bounded, Bounded)) - slope(x, y) = (5 + tanh(40*(x - 1/6)) + tanh(40*(x - 2/6)) + tanh(40*(x - 3/6)) + tanh(40*(x - 4/6)) + tanh(40*(x - 5/6))) / 20 + + slope(x, y) = (5 + tanh(40*(x - 1/6)) + tanh(40*(x - 2/6)) + tanh(40*(x - 3/6)) + tanh(40*(x - 4/6)) + tanh(40*(x - 5/6))) / 20 + (5 + tanh(40*(y - 1/6)) + tanh(40*(y - 2/6)) + tanh(40*(y - 3/6)) + tanh(40*(y - 4/6)) + tanh(40*(y - 5/6))) / 20 # grid = ImmersedBoundaryGrid(grid, GridFittedBottom(slope)) @@ -75,14 +75,14 @@ function setup_simulation(model) Δt = 1e-3 simulation = Simulation(model; Δt = Δt, stop_time = 10, minimum_relative_step = 1e-10) conjure_time_step_wizard!(simulation, cfl=0.7, IterationInterval(1)) - + wall_time = Ref(time_ns()) d = Field{Center, Center, Center}(grid) function progress(sim) pressure_solver = sim.model.pressure_solver - + if pressure_solver isa ConjugateGradientPoissonSolver pressure_iters = iteration(pressure_solver) else @@ -91,11 +91,11 @@ function setup_simulation(model) msg = @sprintf("iter: %d, time: %s, Δt: %.4f, Poisson iters: %d", iteration(sim), prettytime(time(sim)), sim.Δt, pressure_iters) - + elapsed = 1e-9 * (time_ns() - wall_time[]) - + compute_flow_divergence!(d, sim.model) - + msg *= @sprintf(", max u: %6.3e, max v: %6.3e, max w: %6.3e, max b: %6.3e, max d: %6.3e, max pressure: %6.3e, wall time: %s", maximum(sim.model.velocities.u), maximum(sim.model.velocities.v), @@ -104,10 +104,10 @@ function setup_simulation(model) maximum(d), maximum(sim.model.pressures.pNHS), prettytime(elapsed)) - + @info msg wall_time[] = time_ns() - + return nothing end @@ -146,7 +146,7 @@ function setup_simulation(model) filename = filename, schedule = TimeInterval(0.1), overwrite_existing = true) - + return simulation, file_prefix end diff --git a/validation/immersed_boundaries/tracer_advection_over_bump.jl b/validation/immersed_boundaries/tracer_advection_over_bump.jl index 81eb8a980b2..7785684b673 100644 --- a/validation/immersed_boundaries/tracer_advection_over_bump.jl +++ b/validation/immersed_boundaries/tracer_advection_over_bump.jl @@ -105,4 +105,3 @@ simulation.output_writers[:fields] = JLD2Writer(model, model.tracers, force = true) run!(simulation) - diff --git a/validation/immersed_boundaries/vortex_sheet.jl b/validation/immersed_boundaries/vortex_sheet.jl index 472acfa5945..6089a443b35 100644 --- a/validation/immersed_boundaries/vortex_sheet.jl +++ b/validation/immersed_boundaries/vortex_sheet.jl @@ -201,4 +201,4 @@ display(fig) CairoMakie.record(fig, "./vortex_sheet.mp4", 1:Nt, framerate=15, px_per_unit=2) do nn n[] = nn end -# #%% \ No newline at end of file +# #%% diff --git a/validation/implicit_dissipation/variance_dissipation.jl b/validation/implicit_dissipation/variance_dissipation.jl index b540495b6fd..63f94d27194 100644 --- a/validation/implicit_dissipation/variance_dissipation.jl +++ b/validation/implicit_dissipation/variance_dissipation.jl @@ -59,15 +59,15 @@ c⁻ = CenterField(grid) Δtc² = CenterField(grid) for (ts, timestepper) in zip((:AB2, :RK3), (:QuasiAdamsBashforth2, :SplitRungeKutta3)) - - model = HydrostaticFreeSurfaceModel(grid; - timestepper, - velocities, - tracer_advection, - closure, + + model = HydrostaticFreeSurfaceModel(grid; + timestepper, + velocities, + tracer_advection, + closure, tracers=:c, auxiliary_fields=(; Δtc², c⁻)) - + set!(model, c=c₀) set!(model.auxiliary_fields.c⁻, c₀) @@ -91,7 +91,7 @@ for (ts, timestepper) in zip((:AB2, :RK3), (:QuasiAdamsBashforth2, :SplitRungeKu overwrite_existing=true) sim.callbacks[:compute_tracer_dissipation] = Callback(compute_tracer_dissipation!, IterationInterval(1)) - + run!(sim) end @@ -133,4 +133,4 @@ scatter!(ax, rtimes, r_∫closs, label="RK3 total variance loss", color=:blue, m lines!(ax, rtimes, r_∫A, label="RK3 advection dissipation", color=:red, linestyle=:dash) lines!(ax, rtimes, r_∫D, label="RK3 diffusive dissipation", color=:green, linestyle=:dash) lines!(ax, rtimes, r_∫T, label="RK3 total dissipation", color=:purple, linestyle=:dash) -Legend(fig[1, 2], ax) \ No newline at end of file +Legend(fig[1, 2], ax) diff --git a/validation/implicit_free_surface/geostrophic_adjustment_test.jl b/validation/implicit_free_surface/geostrophic_adjustment_test.jl index c9d6b2335e1..b4812d49b57 100644 --- a/validation/implicit_free_surface/geostrophic_adjustment_test.jl +++ b/validation/implicit_free_surface/geostrophic_adjustment_test.jl @@ -176,4 +176,4 @@ function plot_variable2(sims, var1, var2; @info "Frame $i of $Nt" iter[] = i end -end \ No newline at end of file +end diff --git a/validation/implicit_free_surface/nonhydrostatic_implicit_free_surface.jl b/validation/implicit_free_surface/nonhydrostatic_implicit_free_surface.jl index fb14f135ea3..937bee5284d 100644 --- a/validation/implicit_free_surface/nonhydrostatic_implicit_free_surface.jl +++ b/validation/implicit_free_surface/nonhydrostatic_implicit_free_surface.jl @@ -29,7 +29,7 @@ c = sqrt(g) simulation = Simulation(model; Δt, stop_iteration=10) #stop_time=5/c) ηt = [] -function progress(sim) +function progress(sim) @info @sprintf("Time: %s, iteration: %d", prettytime(sim), iteration(sim)) push!(ηt, deepcopy(interior(model.free_surface.displacement, :, 1, 1))) return nothing diff --git a/validation/isopycnal_skew_symmetric_diffusion/coarse_baroclinic_adjustment.jl b/validation/isopycnal_skew_symmetric_diffusion/coarse_baroclinic_adjustment.jl index da54851a7ba..13edf3b98b2 100644 --- a/validation/isopycnal_skew_symmetric_diffusion/coarse_baroclinic_adjustment.jl +++ b/validation/isopycnal_skew_symmetric_diffusion/coarse_baroclinic_adjustment.jl @@ -103,7 +103,7 @@ function progress(sim) prettytime(sim.Δt)) wall_clock[] = time_ns() - + return nothing end @@ -186,4 +186,3 @@ record(fig, filename * ".mp4", 1:Nt, framerate=8) do i @info "Plotting frame $i of $Nt" n[] = i end - diff --git a/validation/isopycnal_skew_symmetric_diffusion/coarse_lat_lon_baroclinic_adjustment.jl b/validation/isopycnal_skew_symmetric_diffusion/coarse_lat_lon_baroclinic_adjustment.jl index 2d5453ef7f9..fed0f3a683a 100644 --- a/validation/isopycnal_skew_symmetric_diffusion/coarse_lat_lon_baroclinic_adjustment.jl +++ b/validation/isopycnal_skew_symmetric_diffusion/coarse_lat_lon_baroclinic_adjustment.jl @@ -191,4 +191,3 @@ record(fig, filename * ".mp4", 1:Nt, framerate=8) do i @info "Plotting frame $i of $Nt" n[] = i end - diff --git a/validation/isopycnal_skew_symmetric_diffusion/new_front_relax.jl b/validation/isopycnal_skew_symmetric_diffusion/new_front_relax.jl index df7a84c8927..5e30d65d479 100644 --- a/validation/isopycnal_skew_symmetric_diffusion/new_front_relax.jl +++ b/validation/isopycnal_skew_symmetric_diffusion/new_front_relax.jl @@ -29,7 +29,7 @@ viscAh=300. viscAz=2.e-4 grid = RectilinearGrid(topology = (Flat, Bounded, Bounded), - size = (Ny, Nz), + size = (Ny, Nz), y = (-Ly/2, Ly/2), z = ([0 cumsum(reverse(dz))']'.-Lz), halo = (3, 3)) @@ -98,7 +98,7 @@ function progress(sim) prettytime(sim.Δt)) wall_clock[] = time_ns() - + return nothing end @@ -175,4 +175,3 @@ record(fig, filename * ".mp4", 1:Nt, framerate=8) do i @info "Plotting frame $i of $Nt" n[] = i end - diff --git a/validation/isopycnal_skew_symmetric_diffusion/zonally_averaged_baroclinic_adjustment.jl b/validation/isopycnal_skew_symmetric_diffusion/zonally_averaged_baroclinic_adjustment.jl index acd68150fa1..310a6e22d69 100644 --- a/validation/isopycnal_skew_symmetric_diffusion/zonally_averaged_baroclinic_adjustment.jl +++ b/validation/isopycnal_skew_symmetric_diffusion/zonally_averaged_baroclinic_adjustment.jl @@ -218,4 +218,3 @@ record(fig, filename * ".mp4", 1:Nt, framerate=8) do i @info "Plotting frame $i of $Nt" n[] = i end - diff --git a/validation/lagrangian_particles/drogued_buoy.jl b/validation/lagrangian_particles/drogued_buoy.jl index 70aaed8fe0b..626285fed86 100644 --- a/validation/lagrangian_particles/drogued_buoy.jl +++ b/validation/lagrangian_particles/drogued_buoy.jl @@ -15,9 +15,9 @@ drogue_depths = -20:20/(n_particles-1):0 c = CenterField(grid) -particle_properties = StructArray{CTrackingParticle}((zeros(n_particles), - zeros(n_particles), - zeros(n_particles), +particle_properties = StructArray{CTrackingParticle}((zeros(n_particles), + zeros(n_particles), + zeros(n_particles), zeros(n_particles))) particles = LagrangianParticles(particle_properties; dynamics = DroguedParticleDynamics(drogue_depths), tracked_fields = (; c)) @@ -56,8 +56,8 @@ simulation = Simulation(model, Δt = 10, stop_time = 4hours) conjure_time_step_wizard!(simulation, IterationInterval(100), cfl = 0.5, diffusive_cfl = 0.5) simulation.output_writers[:velocities] = JLD2Writer(model, model.velocities; - overwrite_existing = true, - filename = "drogued_velocities.jld2", + overwrite_existing = true, + filename = "drogued_velocities.jld2", schedule = TimeInterval(10minutes)) prog(sim) = @info prettytime(sim) * " in " * prettytime(sim.run_wall_time) * " with Δt = " * prettytime(sim.Δt) @@ -67,16 +67,16 @@ add_callback!(simulation, prog, IterationInterval(50)) run!(simulation) simulation.output_writers[:particles] = JLD2Writer(model, (; particles); - overwrite_existing = true, - filename = "drogued_particles.jld2", + overwrite_existing = true, + filename = "drogued_particles.jld2", schedule = TimeInterval(0.1minutes)) simulation.output_writers[:tracer] = JLD2Writer(model, model.tracers; - overwrite_existing = true, - filename = "drogued_tracer.jld2", + overwrite_existing = true, + filename = "drogued_tracer.jld2", schedule = TimeInterval(0.1minutes)) -particles.properties.x .= 0 -particles.properties.y .= 0 +particles.properties.x .= 0 +particles.properties.y .= 0 simulation.stop_time += 0.2hours @@ -110,4 +110,3 @@ heatmap!(ax, xnodes(c), ynodes(c), c_plt, colorrange = (0, 1), alpha = 0.5) scatter!(ax, x_plt, y_plt, color = particle_c_plt, colorrange = (0, 1)) heatmap!(ax2, xnodes(c), znodes(c), c_slice_plt, colorrange = (0, 1), alpha = 0.5) scatter!(ax2, x_plt, drogue_depths, color = particle_c_plt, colorrange = (0, 1)) - diff --git a/validation/large_eddy_simulation/decaying_homogeneous_isotropic_turbulence.jl b/validation/large_eddy_simulation/decaying_homogeneous_isotropic_turbulence.jl index 93cd86ffd0b..01c0af074d1 100644 --- a/validation/large_eddy_simulation/decaying_homogeneous_isotropic_turbulence.jl +++ b/validation/large_eddy_simulation/decaying_homogeneous_isotropic_turbulence.jl @@ -59,11 +59,11 @@ add_callback!(simulation, add_line!, TimeInterval(0.1)) run!(simulation) lines!(ax, [exp(1), exp(3)], x->(x/exp(1))^(-5/3)*exp(15), color = :red, linestyle = :dash) -text!(ax, exp(0.8), exp(13); text = L"$E(k)\sim k^{-5/3}$", color = :white) +text!(ax, exp(0.8), exp(13); text = L"$E(k)\sim k^{-5/3}$", color = :white) Colorbar(fig[1, 2], colormap = :oslo, colorrange = (0, 10), label = "Time (s)") k_filt = 1/sqrt(closure.Cν * 3 / (1/(2*minimum_xspacing(grid))^2 + 1/(2*minimum_yspacing(grid))^2 + 1/(2*minimum_zspacing(grid))^2)) lines!(ax, ones(2) .* k_filt, [exp(0), exp(10)], color = :red, linestyle = :dash) text!(ax, k_filt, exp(10); text = "1/δ") -save("energy.png", fig) \ No newline at end of file +save("energy.png", fig) diff --git a/validation/lid_driven_cavity/lid_driven_cavity.jl b/validation/lid_driven_cavity/lid_driven_cavity.jl index 4ec1c7e06b4..241c61efffd 100644 --- a/validation/lid_driven_cavity/lid_driven_cavity.jl +++ b/validation/lid_driven_cavity/lid_driven_cavity.jl @@ -69,4 +69,3 @@ simulate_lid_driven_cavity(Re=3200, N=128, end_time=50) simulate_lid_driven_cavity(Re=5000, N=256, end_time=50) simulate_lid_driven_cavity(Re=7500, N=256, end_time=75) simulate_lid_driven_cavity(Re=10000, N=256, end_time=100) - diff --git a/validation/mesoscale_turbulence/coarse_function.jl b/validation/mesoscale_turbulence/coarse_function.jl index 2f2c203c7a2..ad0c128f098 100644 --- a/validation/mesoscale_turbulence/coarse_function.jl +++ b/validation/mesoscale_turbulence/coarse_function.jl @@ -154,4 +154,4 @@ for i in eachindex(clist) end end -# (; κskew = 2000.0, κsymmetric = 2000.0, tapering = 0.01, scale = 2.0) failed \ No newline at end of file +# (; κskew = 2000.0, κsymmetric = 2000.0, tapering = 0.01, scale = 2.0) failed diff --git a/validation/multi_region/multi_region_near_global_quarter_degree.jl b/validation/multi_region/multi_region_near_global_quarter_degree.jl index 5dfa16d29ed..21845bb0e4f 100644 --- a/validation/multi_region/multi_region_near_global_quarter_degree.jl +++ b/validation/multi_region/multi_region_near_global_quarter_degree.jl @@ -320,5 +320,3 @@ run!(simulation, pickup = pickup_file) Free surface: $(typeof(model.free_surface).name.wrapper) Time step: $(prettytime(Δt)) """ - - diff --git a/validation/multi_region/multi_region_turbulence.jl b/validation/multi_region/multi_region_turbulence.jl index e8da1929c50..acf1c5c8e1b 100644 --- a/validation/multi_region/multi_region_turbulence.jl +++ b/validation/multi_region/multi_region_turbulence.jl @@ -86,4 +86,3 @@ elapsed_time = 1e-9 * (time_ns() - start_time) u_2, v_2, w = model_2.velocities ζ_2 = compute!(Field(∂x(v_2) - ∂y(u_2))) - diff --git a/validation/open_boundaries/2d_flow_over_bathymetry.jl b/validation/open_boundaries/2d_flow_over_bathymetry.jl index f087b27466e..c52ee152454 100644 --- a/validation/open_boundaries/2d_flow_over_bathymetry.jl +++ b/validation/open_boundaries/2d_flow_over_bathymetry.jl @@ -180,4 +180,4 @@ common_kwargs = (; arch = CPU(), simulation = create_mass_conservation_simulation(; use_open_boundary_condition = true, animation = true, common_kwargs...); run!(simulation) -save("2d_flow_over_bathymetry.mp4", io) \ No newline at end of file +save("2d_flow_over_bathymetry.mp4", io) diff --git a/validation/open_boundaries/oscillating_flow.jl b/validation/open_boundaries/oscillating_flow.jl index d0a43a31732..d9817a6c0cd 100644 --- a/validation/open_boundaries/oscillating_flow.jl +++ b/validation/open_boundaries/oscillating_flow.jl @@ -146,4 +146,3 @@ for grid in (xygrid, xzgrid) run_cylinder(grid, boundary_conditions, simname = simname, stop_time = T) end end - diff --git a/validation/orthogonal_spherical_shell_grid/polar_turbulence.jl b/validation/orthogonal_spherical_shell_grid/polar_turbulence.jl index 103a3e78555..9de53bd4f87 100644 --- a/validation/orthogonal_spherical_shell_grid/polar_turbulence.jl +++ b/validation/orthogonal_spherical_shell_grid/polar_turbulence.jl @@ -105,4 +105,3 @@ frames = 1:length(times) record(fig, "polar_turbulence.mp4", frames, framerate = 12) do i n[] = i end - diff --git a/validation/orthogonal_spherical_shell_grid/rotated_lat_lon_splash.jl b/validation/orthogonal_spherical_shell_grid/rotated_lat_lon_splash.jl index 73fd67917dd..1afc9c09aaa 100644 --- a/validation/orthogonal_spherical_shell_grid/rotated_lat_lon_splash.jl +++ b/validation/orthogonal_spherical_shell_grid/rotated_lat_lon_splash.jl @@ -129,4 +129,3 @@ frames = 1:length(times) record(fig, "splash.mp4", frames, framerate = 12) do nn n[] = nn end - diff --git a/validation/regridding/rectilinear_regridding.jl b/validation/regridding/rectilinear_regridding.jl index 925ae9014ad..d865043770c 100644 --- a/validation/regridding/rectilinear_regridding.jl +++ b/validation/regridding/rectilinear_regridding.jl @@ -55,4 +55,3 @@ scatter!(ax5, Array(interior(c_xy, 180, 60, :)), z) scatter!(ax6, Array(interior(c_xyz, 180, 60, :)), z′) display(fig) - diff --git a/validation/solid_body_rotation/rossby_haurwitz.jl b/validation/solid_body_rotation/rossby_haurwitz.jl index 1e32bfc73ca..8953e8e13ef 100644 --- a/validation/solid_body_rotation/rossby_haurwitz.jl +++ b/validation/solid_body_rotation/rossby_haurwitz.jl @@ -125,4 +125,3 @@ end filepath_w = run_rossby_haurwitz(architecture=GPU(), Nx=512, Ny=256, advection_scheme=WENO(vector_invariant=VelocityStencil()), prefix = "WENOVectorInvariantVel") filepath_w = run_rossby_haurwitz(architecture=GPU(), Nx=512, Ny=256, advection_scheme=WENO(vector_invariant=VorticityStencil()), prefix = "WENOVectorInvariantVort") filepath_w = run_rossby_haurwitz(architecture=GPU(), Nx=512, Ny=256, advection_scheme=VectorInvariant(), prefix = "VectorInvariant") - diff --git a/validation/solid_body_rotation/visualization.jl b/validation/solid_body_rotation/visualization.jl index ef6bf7a35ee..8791b02f165 100644 --- a/validation/solid_body_rotation/visualization.jl +++ b/validation/solid_body_rotation/visualization.jl @@ -159,4 +159,4 @@ function visualize_cartesian_comparison(filepath1, filepath2, title1, title2, va end return nothing -end \ No newline at end of file +end diff --git a/validation/stokes_drift/Langmuir_with_Stokes_x_jet.jl b/validation/stokes_drift/Langmuir_with_Stokes_x_jet.jl index b34a55ea5db..ae16b0db725 100644 --- a/validation/stokes_drift/Langmuir_with_Stokes_x_jet.jl +++ b/validation/stokes_drift/Langmuir_with_Stokes_x_jet.jl @@ -501,4 +501,3 @@ end nothing #hide # ![](Stokes_drift_x_jet.mp4) - diff --git a/validation/stokes_drift/surface_wave_quasi_geostrophic_flow.jl b/validation/stokes_drift/surface_wave_quasi_geostrophic_flow.jl index d5854b22aa3..26299ff9595 100644 --- a/validation/stokes_drift/surface_wave_quasi_geostrophic_flow.jl +++ b/validation/stokes_drift/surface_wave_quasi_geostrophic_flow.jl @@ -161,4 +161,3 @@ contour!(axw, xc, yc, An, color=:gray, levels=5) record(fig, "surface_wave_quasi_geostrophic_induced_flow.mp4", 1:Nt, framerate=8) do nn n[] = nn end - diff --git a/validation/stommel_gyre/stommel_gyre.jl b/validation/stommel_gyre/stommel_gyre.jl index 2a3e9ca1c6d..be936121b07 100644 --- a/validation/stommel_gyre/stommel_gyre.jl +++ b/validation/stommel_gyre/stommel_gyre.jl @@ -14,4 +14,4 @@ integrated velocity. @inline u_Stommel_max(; L, τ₀, β) = τ₀*π^2/β/L @inline v_Stommel_max(; L, τ₀, β) = τ₀*π/β/L -@inline ϵ(r, β, L) = r / (β*L) \ No newline at end of file +@inline ϵ(r, β, L) = r / (β*L) diff --git a/validation/stratified_couette_flow/run_stratified_couette_flow_simulations.jl b/validation/stratified_couette_flow/run_stratified_couette_flow_simulations.jl index 89d90d61236..c1d5fed1ec2 100644 --- a/validation/stratified_couette_flow/run_stratified_couette_flow_simulations.jl +++ b/validation/stratified_couette_flow/run_stratified_couette_flow_simulations.jl @@ -3,4 +3,3 @@ include("stratified_couette_flow.jl") simulate_stratified_couette_flow(Nxy=128, Nz=128, Ri=0) simulate_stratified_couette_flow(Nxy=128, Nz=128, Ri=0.01) simulate_stratified_couette_flow(Nxy=128, Nz=128, Ri=0.04) - diff --git a/validation/stratified_couette_flow/stratified_couette_flow_analysis.jl b/validation/stratified_couette_flow/stratified_couette_flow_analysis.jl index 3bb2ba7bf78..33a89e76381 100644 --- a/validation/stratified_couette_flow/stratified_couette_flow_analysis.jl +++ b/validation/stratified_couette_flow/stratified_couette_flow_analysis.jl @@ -266,4 +266,3 @@ png_filepath = "stratified_couette_flow_velocity_temperature_slices.png" @info "Saving $png_filepath..." savefig(png_filepath, dpi=200) close(fig) - diff --git a/validation/vertical_mixing_closures/bottom_boundary_layer.jl b/validation/vertical_mixing_closures/bottom_boundary_layer.jl index 449831ef075..ab15485383c 100644 --- a/validation/vertical_mixing_closures/bottom_boundary_layer.jl +++ b/validation/vertical_mixing_closures/bottom_boundary_layer.jl @@ -113,4 +113,3 @@ display(fig) record(fig, "bottom_boundary_layer.mp4", 1:Nt, framerate=12) do nn n[] = nn end - diff --git a/validation/vertical_mixing_closures/compare_catke_results.jl b/validation/vertical_mixing_closures/compare_catke_results.jl index 4d0fc3c7e44..d5e5fae25a9 100644 --- a/validation/vertical_mixing_closures/compare_catke_results.jl +++ b/validation/vertical_mixing_closures/compare_catke_results.jl @@ -142,4 +142,3 @@ display(fig) record(fig, "compare_catkes_windy_convection.mp4", 1:Nt, framerate=24) do nn n[] = nn end - diff --git a/validation/vertical_mixing_closures/look_at_stability_functions.jl b/validation/vertical_mixing_closures/look_at_stability_functions.jl index 88f8afe025c..2bee04f9805 100644 --- a/validation/vertical_mixing_closures/look_at_stability_functions.jl +++ b/validation/vertical_mixing_closures/look_at_stability_functions.jl @@ -64,4 +64,3 @@ cf = contourf!(ax6, αᴺ, αᴹ, Pr, levels=0.3:0.1:3.0, colorrrange=(0.35, 2.8 Colorbar(fig[3, 4], cf, vertical=false, tellwidth=false, label="Prandtl number", flipaxis=false) display(fig) - diff --git a/validation/vertical_mixing_closures/test_single_column_model.jl b/validation/vertical_mixing_closures/test_single_column_model.jl index 84a967401da..4c72de300aa 100644 --- a/validation/vertical_mixing_closures/test_single_column_model.jl +++ b/validation/vertical_mixing_closures/test_single_column_model.jl @@ -25,4 +25,3 @@ time_step!(model, 600) @time for n = 1:100 time_step!(model, 600) end - diff --git a/validation/z_star_coordinate/lock_release.jl b/validation/z_star_coordinate/lock_release.jl index a4cbe2a5f9a..75345515b86 100644 --- a/validation/z_star_coordinate/lock_release.jl +++ b/validation/z_star_coordinate/lock_release.jl @@ -21,7 +21,7 @@ model = HydrostaticFreeSurfaceModel(grid; tracers = (:b, :c), timestepper = :SplitRungeKutta3, vertical_coordinate = ZStarCoordinate(grid), - free_surface = SplitExplicitFreeSurface(grid; substeps=20)) # + free_surface = SplitExplicitFreeSurface(grid; substeps=20)) # g = model.free_surface.gravitational_acceleration bᵢ(x, z) = x > 32kilometers ? 0.06 : 0.01 @@ -77,7 +77,7 @@ end simulation.callbacks[:progress] = Callback(progress, IterationInterval(100)) run!(simulation) - + b = model.tracers.b x, y, z = nodes(b) @@ -87,6 +87,6 @@ lines!(ax, (vav .- vav[1]) ./ vav[1], label = "Volume anomaly") lines!(ax, (bav .- bav[1]) ./ bav[1], label = "Buoyancy anomaly") lines!(ax, (cav .- cav[1]) ./ cav[1], label = "Tracer anomaly") axislegend(ax, position=:lt) -ax = Axis(fig[1, 2], title = "Final buoyancy field") +ax = Axis(fig[1, 2], title = "Final buoyancy field") contourf!(ax, x, z, interior(model.tracers.b, :, 1, :), colormap=:balance, levels=20) vlines!(ax, 62.3e3, linestyle = :dash, linewidth = 3, color = :black)