add a unit test for background flux divergence and clean up unrelated commits

[liuchihl]

This PR adds a unit test for #3646

The tests compare model results with and without background fields using a linear stratification profile.

[liuchihl]

@glwagner, this PR has been updated and aligned with the target branch.

[liuchihl]

This unit test has been updated to a more robust correctness test.
Here is a summary of this test:
The test compares two model runs: one with background fields and one without. The model with background fields uses a BackgroundField for buoyancy, while the model without background fields incorporates the stratification as an initial condition. The test then compares the total buoyancy fields (B) from both runs. This verifies that the background field implementation produces consistent results with the traditional approach of including background stratification in the initial conditions.

[glwagner]

Is it necessary to use output here? I think you can directly return the value of B without output. That will make the test simpler and easier to maintain which is super important.

[glwagner]

Main feedback is to set up the test without using an output writer which I don't think is necessary for this test. This will make the test easier to maintain.

Once that is done let me know and I will run this manually and then merge it.

[liuchihl]

Thanks, @glwagner for reviewing this!
I agree with you that it is best to avoid outputwriter in a unit test, but I have trouble reading B:

julia> b = model.tracers.b
1×1×10 Field{Center, Center, Center} on RectilinearGrid on CPU
├── grid: 1×1×10 RectilinearGrid{Float64, Flat, Flat, Bounded} on CPU with 0×0×3 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Nothing, east: Nothing, south: Nothing, north: Nothing, bottom: Gradient, top: Gradient, immersed: ZeroFlux
└── data: 1×1×16 OffsetArray(::Array{Float64, 3}, 1:1, 1:1, -2:13) with eltype Float64 with indices 1:1×1:1×-2:13
    └── max=0.0, min=0.0, mean=0.0

julia> B̄ = model.background_fields.tracers.b
FunctionField located at (Center, Center, Center)
├── func: constant_stratification (generic function with 1 method)
├── grid: 1×1×10 RectilinearGrid{Float64, Flat, Flat, Bounded} on CPU with 0×0×3 halo
├── clock: Clock{Float64, Float64}(time=0 seconds, iteration=0, last_Δt=Inf days)
└── parameters: (N² = 1.0e-6,)

julia> B = B̄ + b 
BinaryOperation at (Center, Center, Center)
├── grid: 1×1×10 RectilinearGrid{Float64, Flat, Flat, Bounded} on CPU with 0×0×3 halo
└── tree: 
    + at (Center, Center, Center)
    ├── 1×1×10 Oceananigans.Fields.FunctionField{Center, Center, Center, Clock{Float64, Float64}, @NamedTuple{N²::Float64}, typeof(constant_stratification), RectilinearGrid{Float64, Flat, Flat, Bounded, Oceananigans.Grids.StaticVerticalCoordinate{OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, Float64}, Float64, Float64, Nothing, Nothing, CPU}, Float64}
    └── 1×1×10 Field{Center, Center, Center} on RectilinearGrid on CPU

As shown, B becomes a BinaryOperation when combining the background and perturbation buoyancy. I have trouble reading B properly without using outputwriter. This contrasts with Field, where I can simply utilize the interior function. Do you have suggestions on how to do this?