Tracer Conservation test for `TripolarGrid` by taimoorsohail · Pull Request #346 · NumericalEarth/NumericalEarth.jl

taimoorsohail · 2026-06-12T02:20:32Z

This PR adds a tracer conservation test for a TripolarGrid setup. Note that at the moment the test is quite heavy. I would appreciate your advice in making it a bit less heavy, but I wasn't sure what components are important (or not) in ensuring conservation. We know tracers are conserved in Oceananigans.jl, so if there is an issue it is likely in the coupling with other datasets. Hence I took a maximalist approach. Let's work on making this PR more lightweight but still revealing if there is a budget non-closure.

codecov · 2026-06-12T02:32:10Z

Codecov Report

✅ All modified and coverable lines are covered by tests.

📢 Thoughts on this report? Let us know!

taimoorsohail · 2026-06-12T03:04:45Z

The checks are:

Heat / OHC

$$ \mathrm{OHC}_n - \mathrm{OHC}_{n-1} \approx -\left( \int_A Q_{n-1} dV \right) (t_n - t_{n-1}) $$

where

$$\mathrm{OHC} = \rho_0 c_p \int_V T dV$$

and $Q$ is the diagnosed net ocean heat flux.

Freshwater content

$$\mathrm{FW}_n - \mathrm{FW}_{n-1} \approx -\left( \int_A F_{n-1} dV \right) (t_n - t_{n-1}) $$

where

$$\mathrm{FW} = -\frac{\rho_0}{S_0} \int_V S dV$$

and $F$ is the diagnosed net ocean freshwater flux.

Sign convention

The minus sign is required because the diagnosed surface fluxes are defined with the opposite sign to the content tendency:

negative heat flux into the ocean corresponds to increasing OHC
negative freshwater flux into the ocean corresponds to increasing freshwater content

Pass criterion

For each timestep, the test requires:

ΔOHC ≈ - previous_heat_flux * Δt
ΔFW ≈ - previous_fw_flux * Δt

to within numerical tolerance.

…lEarth/NumericalEarth.jl into ts/add-tracer-conservation-test

navidcy · 2026-06-12T04:49:20Z

+            ohc_integral = Field(Integral(simulation.model.ocean.model.tracers.T, dims=(1, 2, 3)))
+            fw_content_integral = Field(Integral(simulation.model.ocean.model.tracers.S, dims=(1, 2, 3)))


Suggested change

ohc_integral = Field(Integral(simulation.model.ocean.model.tracers.T, dims=(1, 2, 3)))

fw_content_integral = Field(Integral(simulation.model.ocean.model.tracers.S, dims=(1, 2, 3)))

ocean_heat_content = Field(ρᵒᶜ * cᵒᶜ * Integral(simulation.model.ocean.model.tracers.T, dims=(1, 2, 3)))

freshwater_content = Field(-ρᵒᶜ / Sᵒᶜ * Integral(simulation.model.ocean.model.tracers.S, dims=(1, 2, 3)))

navidcy · 2026-06-12T04:49:48Z

+            heat_flux_integral = Field(Integral(surface_forcing.heat_flux, dims=(1, 2, 3)))
+            fw_flux_integral = Field(Integral(surface_forcing.fw_flux, dims=(1, 2, 3)))


Why are these 3D integrals? Aren't these 2D fields?

Also, when we integrate the heat flux it becomes heat rate, right?

Suggested change

heat_flux_integral = Field(Integral(surface_forcing.heat_flux, dims=(1, 2, 3)))

fw_flux_integral = Field(Integral(surface_forcing.fw_flux, dims=(1, 2, 3)))

heat_rate = Field(Integral(surface_forcing.heat_flux, dims=(1, 2, 3)))

freshwater_rate = Field(Integral(surface_forcing.fw_flux, dims=(1, 2, 3)))

navidcy · 2026-06-12T04:54:59Z

+                current_ohc = ρᵒᶜ * cᵒᶜ * scalar_integral(ohc_integral)
+                current_fw_content = -ρᵒᶜ / Sᵒᶜ * scalar_integral(fw_content_integral)
+                current_time = Float64(simulation.model.clock.time)
+                Δt_model = current_time - previous_time


isn't this the Δt? Why define new variable?
Perhaps just a @test current_time - previous_time == Δt?

To make it clear that dt = t_n - t_n-1 rather than dt = t_n+1 - t_n

What's the difference? The simulation you set up has a fixed Δt anyway?

perhaps we need to call it previous_Δt?

It's fixed here but perhaps not in all cases... yes previous_dt makes sense

I think previous_Δt is saved in clock... So you can use model.clock.last_Δt?

Updated grid size and modified initialization of ocean and sea ice models to include date parameter.

simone-silvestri · 2026-06-12T06:10:38Z

For the conservation to be machine precision on a RightCenterFolded tripolar grid we need CliMA/Oceananigans.jl#5099 which should be ready to merge (it has been hanging for a little bit but it should be ready)

glwagner · 2026-06-12T08:00:43Z

+                @test current_ocean_heat_content - previous_ocean_heat_content ≈ -previous_heat_flux * last_Δt rtol=1e-4 atol=1e-2
+                @test current_freshwater_content - previous_freshwater_content ≈ -previous_freshwater_flux * last_Δt rtol=1e-4 atol=1e-2


consider adding a comment to explain these choice of tolerances

definitely
ideally this should pass with machine precision

in the process of debugging we found that even just a HydrostaticFreeSurfaceModel the tracer doesn't change based on the surface flux; I'm preparing a MWE+test for Oceananigans.

hm... I take it back, the test passes in Oceananigans:

include("dependencies_for_runtests.jl") arch = CPU() grid = RectilinearGrid(arch, size=(10, 10, 10), x=(-10, 10), y=(-4, 4), z=(-2, 0), topology = (Bounded, Bounded, Bounded)) Jᶜ_func(x, y, t) = -2.1 + 20 * exp(-x^2 - y^2) Jᶜ = Field{Center, Center, Nothing}(grid) set!(Jᶜ, (x, y) -> Jᶜ_func(x, y, 0)) c_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Jᶜ_func)) model = HydrostaticFreeSurfaceModel(grid, timestepper=:SplitRungeKutta3, tracers=:c, boundary_conditions=(c=c_bcs,), free_surface=SplitExplicitFreeSurface(substeps=15)) set!(model, c = (x, y, z) -> 10 * z + exp(x + y + z)) c = model.tracers.c total_c = Integral(c, dims=(1, 2, 3)) surface_c_flux = Integral(Jᶜ, dims=(1, 2)) Δt = 0.3 simulation = Simulation(model; Δt, stop_iteration=10) for _ = 1:simulation.stop_iteration previous_total_c = first(Field(total_c)) previous_surface_c_flux = first(Field(surface_c_flux)) time_step!(simulation, Δt) current_total_c = first(Field(total_c)) current_surface_c_flux = first(Field(surface_c_flux)) current_time = Float64(simulation.model.clock.time) last_Δt = simulation.model.clock.last_Δt @info "Iteration $(simulation.model.clock.iteration), time = $(current_time) seconds" @show current_total_c @show previous_total_c @show current_total_c - previous_total_c @show - previous_surface_c_flux * last_Δt @test current_total_c - previous_total_c ≈ - previous_surface_c_flux * last_Δt end

which gives

julia> include("test_tracer_conservation.jl") [2026/06/12 16:53:16.882] INFO Initializing simulation... [2026/06/12 16:53:16.883] INFO ... simulation initialization complete (493.791 μs) [2026/06/12 16:53:16.883] INFO Executing initial time step... [2026/06/12 16:53:16.927] INFO ... initial time step complete (43.986 ms). [2026/06/12 16:53:16.939] INFO Iteration 1, time = 0.3 seconds current_total_c = 856826.6742445715 previous_total_c = 856741.5286719211 current_total_c - previous_total_c = 85.14557265036274 -previous_surface_c_flux * last_Δt = 85.14557265106653 [2026/06/12 16:53:17.017] INFO Iteration 2, time = 0.6 seconds current_total_c = 856911.8198172228 previous_total_c = 856826.6742445715 current_total_c - previous_total_c = 85.14557265129406 -previous_surface_c_flux * last_Δt = 85.14557265106653 [2026/06/12 16:53:17.020] INFO Iteration 3, time = 0.8999999999999999 seconds current_total_c = 856996.965389874 previous_total_c = 856911.8198172228 current_total_c - previous_total_c = 85.14557265117764 -previous_surface_c_flux * last_Δt = 85.14557265106653 [2026/06/12 16:53:17.022] INFO Iteration 4, time = 1.2 seconds current_total_c = 857082.1109625254 previous_total_c = 856996.965389874 current_total_c - previous_total_c = 85.14557265141048 -previous_surface_c_flux * last_Δt = 85.14557265106653 [2026/06/12 16:53:17.025] INFO Iteration 5, time = 1.5 seconds current_total_c = 857167.2565351762 previous_total_c = 857082.1109625254 current_total_c - previous_total_c = 85.1455726508284 -previous_surface_c_flux * last_Δt = 85.14557265106653 [2026/06/12 16:53:17.037] INFO Iteration 6, time = 1.8 seconds current_total_c = 857252.4021078274 previous_total_c = 857167.2565351762 current_total_c - previous_total_c = 85.14557265117764 -previous_surface_c_flux * last_Δt = 85.14557265106653 [2026/06/12 16:53:17.040] INFO Iteration 7, time = 2.1 seconds current_total_c = 857337.5476804782 previous_total_c = 857252.4021078274 current_total_c - previous_total_c = 85.1455726508284 -previous_surface_c_flux * last_Δt = 85.14557265106653 [2026/06/12 16:53:17.042] INFO Iteration 8, time = 2.4 seconds current_total_c = 857422.6932531294 previous_total_c = 857337.5476804782 current_total_c - previous_total_c = 85.14557265117764 -previous_surface_c_flux * last_Δt = 85.14557265106653 [2026/06/12 16:53:17.045] INFO Iteration 9, time = 2.6999999999999997 seconds current_total_c = 857507.8388257805 previous_total_c = 857422.6932531294 current_total_c - previous_total_c = 85.14557265117764 -previous_surface_c_flux * last_Δt = 85.14557265106653 [2026/06/12 16:53:17.053] INFO Simulation is stopping after running for 71.352 ms. [2026/06/12 16:53:17.054] INFO Model iteration 10 equals or exceeds stop iteration 10. [2026/06/12 16:53:17.054] INFO Iteration 10, time = 2.9999999999999996 seconds current_total_c = 857592.9843984314 previous_total_c = 857507.8388257805 current_total_c - previous_total_c = 85.1455726508284 -previous_surface_c_flux * last_Δt = 85.14557265106653

Taimoor Sohail and others added 2 commits June 12, 2026 12:17

Added the tracer conservation test

117cf00

Merge branch 'main' into ts/add-tracer-conservation-test

e5b32ac

taimoorsohail requested review from glwagner, navidcy and simone-silvestri June 12, 2026 02:20

Restrict test runs to tracer conservation

95081e5

Taimoor Sohail and others added 5 commits June 12, 2026 13:11

Increase halo size

ca04997

typo

ded46bb

Float32-64

219cffc

Reduce inteprolation steps

98711db

Refactor tracer conservation test for clarity

1060b71

navidcy reviewed Jun 12, 2026

View reviewed changes

Comment thread test/test_tracer_conservation.jl

navidcy reviewed Jun 12, 2026

View reviewed changes

Comment thread test/test_tracer_conservation.jl Outdated

Taimoor Sohail and others added 3 commits June 12, 2026 14:29

import Units

3d93236

Merge branch 'ts/add-tracer-conservation-test' of github.com:Numerica…

5de5e4e

…lEarth/NumericalEarth.jl into ts/add-tracer-conservation-test

Merge branch 'main' into ts/add-tracer-conservation-test

7d931e3

navidcy reviewed Jun 12, 2026

View reviewed changes

Comment thread test/test_tracer_conservation.jl Outdated

navidcy reviewed Jun 12, 2026

View reviewed changes

Comment thread test/test_tracer_conservation.jl Outdated

navidcy reviewed Jun 12, 2026

View reviewed changes

navidcy added 7 commits June 12, 2026 13:16

Adjust grid size and model initialization in tests

49537b7

Updated grid size and modified initialization of ocean and sea ice models to include date parameter.

Update test_tracer_conservation.jl

c260af4

Update test_tracer_conservation.jl

8b35704

cleaner

a81d366

merge main

a92de31

debug

50cc511

fixes

580b141

try to debug

693f026

navidcy reviewed Jun 12, 2026

View reviewed changes

Comment thread test/test_tracer_conservation.jl Outdated

test works, but fails

10f0c47

glwagner reviewed Jun 12, 2026

View reviewed changes

add tracer conservation test

7cbbcc7

		ohc_integral = Field(Integral(simulation.model.ocean.model.tracers.T, dims=(1, 2, 3)))
		fw_content_integral = Field(Integral(simulation.model.ocean.model.tracers.S, dims=(1, 2, 3)))

		heat_flux_integral = Field(Integral(surface_forcing.heat_flux, dims=(1, 2, 3)))
		fw_flux_integral = Field(Integral(surface_forcing.fw_flux, dims=(1, 2, 3)))

		@test current_ocean_heat_content - previous_ocean_heat_content ≈ -previous_heat_flux * last_Δt rtol=1e-4 atol=1e-2
		@test current_freshwater_content - previous_freshwater_content ≈ -previous_freshwater_flux * last_Δt rtol=1e-4 atol=1e-2

Conversation

taimoorsohail commented Jun 12, 2026

Uh oh!

codecov Bot commented Jun 12, 2026

Codecov Report

Uh oh!

taimoorsohail commented Jun 12, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Heat / OHC

Freshwater content

Sign convention

Pass criterion

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

navidcy Jun 12, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

navidcy Jun 12, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

simone-silvestri commented Jun 12, 2026

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

navidcy Jun 12, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants

taimoorsohail commented Jun 12, 2026 •

edited

Loading

navidcy Jun 12, 2026 •

edited

Loading

navidcy Jun 12, 2026 •

edited

Loading

navidcy Jun 12, 2026 •

edited

Loading