Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add support for online remapping #2179

Merged
merged 23 commits into from
Dec 5, 2023
Merged
Changes from 1 commit
Commits
Show all changes
23 commits
Select commit Hold shift + click to select a range
6567c41
Remove outdated comments
Sbozzolo Sep 28, 2023
0f6b56d
Require output_writer in diagnostics
Sbozzolo Nov 14, 2023
0ed3c02
Rework diagnostic writers
Sbozzolo Nov 14, 2023
6f7f7c5
Make time_index first coordinate in NetCDF
Sbozzolo Oct 11, 2023
fb40d8c
Ignore diagnostics with output_every < 0
Sbozzolo Oct 5, 2023
d35e4a6
Add orography to diagnostics
Sbozzolo Nov 20, 2023
1ea3d37
Add elevation profile in NetCDF output
Sbozzolo Oct 23, 2023
92aa93f
Add support for interpolation of topography
Sbozzolo Oct 11, 2023
8cc1a6b
Output default diagnostics by default
Sbozzolo Oct 20, 2023
fadbc2c
Choose timestep that divides 1 day
Sbozzolo Nov 14, 2023
e6c5a84
Help inference of orchestrate_diagnostics
Sbozzolo Oct 30, 2023
0fb1911
Fix NetCDF performance issue (+ more bugs)
Sbozzolo Oct 31, 2023
258a7d4
Remove redunant reduction in output_long_name
Sbozzolo Nov 6, 2023
f157d32
Be more careful with floating point comparisons
Sbozzolo Nov 22, 2023
ec1b50a
Avoid allocations in diagnostics with U/V/WVectors
Sbozzolo Dec 1, 2023
1da73bb
Remove FT from radiation diagnostics
Sbozzolo Nov 30, 2023
73930d7
Bump allocations
Sbozzolo Dec 4, 2023
ab74587
Output relevant diagnostics
Sbozzolo Nov 29, 2023
d69bfe9
Fix YAML parsing of diagnostics:output_name
Sbozzolo Nov 29, 2023
a8bb43b
Add AtmosSimulation
Sbozzolo Nov 29, 2023
27b0932
Fail test with fewer than 0.5 * alloc_limit
Sbozzolo Nov 30, 2023
1558d75
Avoid allocations due to surface_ct3_unit
Sbozzolo Dec 1, 2023
91fe85c
Add warn_allocations_diagnostics
Sbozzolo Dec 1, 2023
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Prev Previous commit
Next Next commit
Fix NetCDF performance issue (+ more bugs)
Sbozzolo committed Dec 4, 2023

Verified

This commit was created on GitHub.com and signed with GitHub’s verified signature.
commit 0fb19114d4281a41a7ed2db5a7d8b10376d5d409
61 changes: 37 additions & 24 deletions src/diagnostics/netcdf_writer.jl
Original file line number Diff line number Diff line change
@@ -69,9 +69,12 @@ function add_dimension_maybe!(
dim.attrib[String(k)] = v
end

# We use copyto! so that we don't have to specify the dimensions (e.g., instead of
# dim[:]), since the dimensions could be different when working with topography.
copyto!(dim, points)
if length(size(points)) == 1
dim[:] = points
else
# We have topography
dim[:, :, :] = points
end
end
return nothing
end
@@ -154,11 +157,6 @@ function add_space_coordinates_maybe!(
return [name]
end

# For the horizontal space, we also have to look at the domain, so we define another set of
# functions that dispatches over the domain
target_coordinates(space::Spaces.AbstractSpectralElementSpace, num_points) =
target_coordinates(space, num_points, Meshes.domain(space.topology.mesh))

add_space_coordinates_maybe!(
nc::NCDatasets.NCDataset,
space::Spaces.AbstractSpectralElementSpace,
@@ -167,9 +165,15 @@ add_space_coordinates_maybe!(
nc,
space,
num_points,
Meshes.domain(space.topology.mesh),
Meshes.domain(space.topology),
)


# For the horizontal space, we also have to look at the domain, so we define another set of
# functions that dispatches over the domain
target_coordinates(space::Spaces.AbstractSpectralElementSpace, num_points) =
target_coordinates(space, num_points, Meshes.domain(space.topology))

# Box
function target_coordinates(
space::Spaces.SpectralElementSpace2D,
@@ -386,16 +390,12 @@ function hcoords_from_horizontal_space(
end

"""
hcoords_from_horizontal_space(space, hpts)
hcoords_from_horizontal_space(space, domain, hpts)

Prepare the matrix of horizontal coordinates with the correct type according to the given
`space` (e.g., `ClimaCore.Geometry.LatLongPoint`s).
Prepare the matrix of horizontal coordinates with the correct type according to the given `space`
and `domain` (e.g., `ClimaCore.Geometry.LatLongPoint`s).
"""
hcoords_from_horizontal_space(space, hpts) = hcoords_from_horizontal_space(
space,
Meshes.domain(space.topology.mesh),
hpts,
)
function hcoords_from_horizontal_space(space, domain, hpts) end

struct NetCDFWriter{T, TS}

@@ -478,7 +478,11 @@ function NetCDFWriter(;
elseif hypsography isa Hypsography.LinearAdaption
horizontal_space = axes(hypsography.surface)
hpts = target_coordinates(horizontal_space, num_points)
hcoords = hcoords_from_horizontal_space(horizontal_space, hpts)
hcoords = hcoords_from_horizontal_space(
horizontal_space,
horizontal_space.topology.mesh.domain,
hpts,
)
vcoords = []
remapper = Remapper(hcoords, vcoords, horizontal_space)
interpolated_surface =
@@ -488,11 +492,11 @@ function NetCDFWriter(;
end

return NetCDFWriter{typeof(num_points), typeof(interpolated_surface)}(
Dict(), # remappers
Dict(),
num_points,
compression_level,
interpolated_surface,
Dict(), # open_files
Dict(),
interpolate_z_over_msl,
)
end
@@ -532,7 +536,11 @@ function write_field!(writer::NetCDFWriter, field, diagnostic, integrator)

# zcoords is going to be empty for a 2D horizontal slice
zcoords = [Geometry.ZPoint(p) for p in vpts]
hcoords = hcoords_from_horizontal_space(horizontal_space, hpts)
hcoords = hcoords_from_horizontal_space(
horizontal_space,
Meshes.domain(horizontal_space.topology),
hpts,
)

writer.remappers[var.short_name] = Remapper(hcoords, zcoords, space)
end
@@ -598,7 +606,12 @@ function write_field!(writer::NetCDFWriter, field, diagnostic, integrator)

nc["time"][time_index] = integrator.t

# selectdim(v, 1, time_index) is equivalent to v[time_index, :, :] or
# v[time_index, :, :, :] in 3/4 dimensions
selectdim(v, 1, time_index) .= interpolated_field
# TODO: It would be nice to find a cleaner way to do this
if length(dim_names) == 3
v[time_index, :, :, :] = interpolated_field
elseif length(dim_names) == 2
v[time_index, :, :] = interpolated_field
elseif length(dim_names) == 1
v[time_index, :] = interpolated_field
end
end