From 89fe52590f2b5e6f82a58b53bc841e249d7fae48 Mon Sep 17 00:00:00 2001 From: Gabriele Bozzola Date: Mon, 2 Oct 2023 13:00:04 -0700 Subject: [PATCH] Rework diagnostic writers MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - Transform writers to standalone objects instead of functions - Have one writer for all the diagnostics - Add support to general NetCDF writer Now with support for all the configurations (including planes ✈️) --- docs/src/diagnostics.md | 41 ++- src/diagnostics/Diagnostics.jl | 7 +- src/diagnostics/diagnostic.jl | 7 +- src/diagnostics/writers.jl | 557 ++++++++++++++++++++++++++++----- src/solver/solve.jl | 3 + src/solver/type_getters.jl | 7 +- 6 files changed, 527 insertions(+), 95 deletions(-) diff --git a/docs/src/diagnostics.md b/docs/src/diagnostics.md index b55ad74b690..8357a66fcfd 100644 --- a/docs/src/diagnostics.md +++ b/docs/src/diagnostics.md @@ -28,9 +28,9 @@ identifies the period over which to compute the reduction and how often to save to disk. `output_name` is optional, and if provided, it identifies the name of the output file. -The default `writer` is HDF5. If `writer` is `nc` or `netcdf`, the output is +The default `writer` is NetCDF. If `writer` is `nc` or `netcdf`, the output is remapped non-conservatively on a Cartesian grid and saved to a NetCDF file. -Currently, only 3D fields on cubed spheres are supported. +Currently, the NetCDF writer does not support plane configurations. ### From a script @@ -91,7 +91,7 @@ More specifically, a `ScheduledDiagnostic` contains the following pieces of data - `variable`: The diagnostic variable that has to be computed and output. - `output_every`: Frequency of how often to save the results to disk. -- `output_writer`: Function that controls out to save the computed diagnostic +- `output_writer`: Object that controls out to save the computed diagnostic variable to disk. - `reduction_time_func`: If not `nothing`, the `ScheduledDiagnostic` receives an area of scratch space `acc` where to accumulate partial results. Then, at @@ -134,18 +134,29 @@ that we have a neutral state for the accumulators that are used in the reductions. A key entry in a `ScheduledDiagnostic` is the `output_writer`, the function -responsible for saving the output to disk. `output_writer` is called with three -arguments: the value that has to be output, the `ScheduledDiagnostic`, and the -integrator. Internally, the integrator contains extra information (such as the -current timestep), so the `output_writer` can implement arbitrarily complex -behaviors. It is responsibility of the `output_writer` to properly use the -provided information for meaningful output. `ClimaAtmos` provides functions that -return `output_writer`s for standard operations. The main one is currently -`HDF5Writer`, which should be enough for most use cases. To use it, just -initialize a `ClimaAtmos.HDF5Writer` object with your choice of configuration -and pass it as a `output_writer` argument to the `ScheduledDiagnostic`. More -information about the options supported by `ClimaAtmos.HDF5Writer` is available -in its constructor. +responsible for saving the output to disk. `output_writer` is an object that +contains information on how to write the diagnostic to disk. It is called by +calling a function `write_field!(writer, field, diagnostic, integrator)`, where +`field` is the value that has to be output, `diagnostic`, the +`ScheduledDiagnostic`, and `integrator` is the `SciML` integrator. Internally, +the integrator contains extra information (such as the current timestep), so the +`output_writer` can implement arbitrarily complex behaviors. It is +responsibility of the `output_writer` to properly use the provided information +for meaningful output. + +`ClimaAtmos` provides `output_writer`s for standard operations. The main writer +is `NetCDFWriter`, which produces NetCDF files for the fields on lat-long-z or +x-y-z coordinates. The `NetCDFWriter` performs a Lagrange interpolation of the +fields onto the desired coordinates. To use a `NetCDFWriter`, you just need to +provide how many points you want the target space to be along the various +dimensions. `NetCDFWriter` will interpolate on a linearly spaced grid defined by +the boundaries of the domain and the number of points provided. + +Another writer one is the `HDF5Writer`, which saves diagnostics as defined on +the computational grid. To use it, just initialize a `ClimaAtmos.HDF5Writer` +object with your choice of configuration and pass it as a `output_writer` +argument to the `ScheduledDiagnostic`. More information about the options +supported by `ClimaAtmos.HDF5Writer` is available in its constructor. There are two flavors of `ScheduledDiagnostic`s: `ScheduledDiagnosticIterations`, and `ScheduledDiagnosticTime`. The main diff --git a/src/diagnostics/Diagnostics.jl b/src/diagnostics/Diagnostics.jl index 17f3d4ff5f7..be3a03ed27a 100644 --- a/src/diagnostics/Diagnostics.jl +++ b/src/diagnostics/Diagnostics.jl @@ -1,8 +1,11 @@ module Diagnostics +import ClimaComms + import LinearAlgebra: dot -import ClimaCore: Fields, Geometry, InputOutput, Meshes, Spaces, Operators +import ClimaCore: + Fields, Geometry, InputOutput, Meshes, Spaces, Operators, Domains import ClimaCore.Utilities: half import Thermodynamics as TD @@ -37,4 +40,6 @@ include(joinpath("..", "utils", "abbreviations.jl")) include("diagnostic.jl") include("writers.jl") + +export close end diff --git a/src/diagnostics/diagnostic.jl b/src/diagnostics/diagnostic.jl index cd83ee7f7d9..57606bf70ee 100644 --- a/src/diagnostics/diagnostic.jl +++ b/src/diagnostics/diagnostic.jl @@ -714,7 +714,12 @@ function get_callbacks_from_diagnostics( diag.pre_output_hook!(storage[diag], counters[diag]) # Write to disk - diag.output_writer(storage[diag], diag, integrator) + write_field!( + diag.output_writer, + storage[diag], + diag, + integrator, + ) # accumulator[diag] is not defined for non-reductions diag_accumulator = get(accumulators, diag, nothing) diff --git a/src/diagnostics/writers.jl b/src/diagnostics/writers.jl index 994a06067a9..2971dfda6c7 100644 --- a/src/diagnostics/writers.jl +++ b/src/diagnostics/writers.jl @@ -3,9 +3,14 @@ # This file defines function-generating functions for output_writers for diagnostics. The # writers come with opinionated defaults. -import ClimaCore.Remapping: interpolate_array +import ClimaCore.Remapping: Remapper, interpolate, interpolate_array + import NCDatasets +############## +# HDF5Writer # +############## + """ HDF5Writer() @@ -29,48 +34,399 @@ We need to implement the following features/options: - ...more features/options """ -function HDF5Writer() - # output_drivers are called with the three arguments: the value, the ScheduledDiagnostic, - # and the integrator - function write_to_hdf5(value, diagnostic, integrator) - var = diagnostic.variable - time = integrator.t +struct HDF5Writer end - output_path = joinpath( - integrator.p.simulation.output_dir, - "$(diagnostic.output_short_name)_$(time).h5", - ) - hdfwriter = InputOutput.HDF5Writer(output_path, integrator.p.comms_ctx) - InputOutput.write!(hdfwriter, value, "$(diagnostic.output_short_name)") - attributes = Dict( - "time" => time, - "long_name" => diagnostic.output_long_name, - "variable_units" => var.units, - "standard_variable_name" => var.standard_name, - ) +""" + close(writer::HDF5Writer) + +Close all the files open in `writer`. (Currently no-op.) +""" +close(writer::HDF5Writer) = nothing + +function write_field!(writer::HDF5Writer, field, diagnostic, integrator) + var = diagnostic.variable + time = integrator.t + + output_path = joinpath( + integrator.p.simulation.output_dir, + "$(diagnostic.output_short_name)_$(time).h5", + ) + + hdfwriter = InputOutput.HDF5Writer(output_path, integrator.p.comms_ctx) + InputOutput.write!(hdfwriter, field, "$(diagnostic.output_short_name)") + attributes = Dict( + "time" => time, + "long_name" => diagnostic.output_long_name, + "variable_units" => var.units, + "standard_variable_name" => var.standard_name, + ) + + # TODO: Use directly InputOutput functions + InputOutput.HDF5.h5writeattr( + hdfwriter.file.filename, + "fields/$(diagnostic.output_short_name)", + attributes, + ) + + Base.close(hdfwriter) + return nothing +end + +################ +# NetCDFWriter # +################ +""" + add_dimension_maybe!(nc::NCDatasets.NCDataset, + name::String, + points::Vector{FT}; + kwargs...) - # TODO: Use directly InputOutput functions - InputOutput.HDF5.h5writeattr( - hdfwriter.file.filename, - "fields/$(diagnostic.output_short_name)", - attributes, - ) - Base.close(hdfwriter) - return nothing +Add dimension identified by `name` in the given `nc` file and fill it with the given +`points`. If the dimension already exists, check if it is consistent with the new one. +Optionally, add all the keyword arguments as attributes. +""" +function add_dimension_maybe!( + nc::NCDatasets.NCDataset, + name::String, + points::Vector{FT}; + kwargs..., +) where {FT <: Real} + + if haskey(nc, "$name") + # dimension already exists: check correct size + if size(nc["$name"]) != size(points) + error("Incompatible $name dimension already exists") + end + else + NCDatasets.defDim(nc, "$name", length(points)) + dim = NCDatasets.defVar(nc, "$name", FT, ("$name",)) + for (k, v) in kwargs + dim.attrib[String(k)] = v + end + dim[:] = points end + return nothing +end + +""" + add_time_maybe!(nc::NCDatasets.NCDataset, + float_type::DataType; + kwargs...) + + +Add the `time` dimension (with infinite size) to the given NetCDF file if not already there. +Optionally, add all the keyword arguments as attributes. +""" +function add_time_maybe!( + nc::NCDatasets.NCDataset, + float_type::DataType; + kwargs..., +) + + # If we already have time, do nothing + haskey(nc, "time") && return nothing + + NCDatasets.defDim(nc, "time", Inf) + dim = NCDatasets.defVar(nc, "time", float_type, ("time",)) + for (k, v) in kwargs + dim.attrib[String(k)] = v + end + return nothing +end + +""" + add_space_coordinates_maybe!(nc::NCDatasets.NCDataset, + space::Spaces.AbstractSpace, + num_points) + +Add dimensions relevant to the `space` to the given `nc` NetCDF file. The range is +automatically determined and the number of points is set with `num_points`, which has to be +an iterable of size N, where N is the number of dimensions of the space. For instance, 3 for +a cubed sphere, 2 for a surface, 1 for a column. + +The function returns an array with the names of the relevant dimensions. (We want arrays +because we want to preserve the order to match the one in num_points). +""" +function add_space_coordinates_maybe! end + +""" + target_coordinates!(space::Spaces.AbstractSpace, + num_points) + +Return the range of interpolation coordinates. The range is automatically determined and the +number of points is set with `num_points`, which has to be an iterable of size N, where N is +the number of dimensions of the space. For instance, 3 for a cubed sphere, 2 for a surface, +1 for a column. +""" +function target_coordinates(space, num_points) end + +function target_coordinates( + space::S, + num_points, +) where { + S <: + Union{Spaces.CenterFiniteDifferenceSpace, Spaces.FaceFiniteDifferenceSpace}, +} + num_points_z = num_points[] + FT = Spaces.undertype(space) + vert_domain = space.topology.mesh.domain + z_min, z_max = vert_domain.coord_min.z, vert_domain.coord_max.z + return collect(range(FT(z_min), FT(z_max), num_points_z)) +end + +# Cell-centered column +function add_space_coordinates_maybe!( + nc::NCDatasets.NCDataset, + space::Spaces.CenterFiniteDifferenceSpace, + num_points_z, +) + name = "z" + zpts = target_coordinates(space, num_points_z) + add_dimension_maybe!(nc, "z", zpts, units = "m") + return [name] +end + +# Face-centered column +function add_space_coordinates_maybe!( + nc::NCDatasets.NCDataset, + space::Spaces.FaceFiniteDifferenceSpace, + num_points_z_half, +) + name = "z_half" + z_half_pts = target_coordinates(space, num_points_z_half) + add_dimension_maybe!(nc, "z_half", z_half_pts, units = "m") + 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, space.topology.mesh.domain) +add_space_coordinates_maybe!( + nc::NCDatasets.NCDataset, + space::Spaces.AbstractSpectralElementSpace, + num_points, +) = add_space_coordinates_maybe!( + nc, + space, + num_points, + space.topology.mesh.domain, +) + +# Box +function target_coordinates( + space::Spaces.SpectralElementSpace2D, + num_points, + domain::Domains.RectangleDomain, +) + num_points_x, num_points_y = num_points + FT = Spaces.undertype(space) + xmin = FT(domain.coord_min.x) + xmax = FT(domain.coord_max.x) + ymin = FT(domain.coord_min.y) + ymax = FT(domain.coord_max.y) + xpts = collect(range(xmin, xmax, num_points_x)) + ypts = collect(range(ymin, ymax, num_points_y)) + + return (xpts, ypts) +end + +# Plane +function target_coordinates( + space::Spaces.SpectralElementSpace1D, + num_points, + domain::Domains.IntervalDomain, +) + num_points_x, _... = num_points + FT = Spaces.undertype(space) + xmin = FT(domain.coord_min.x) + xmax = FT(domain.coord_max.x) + xpts = collect(range(xmin, xmax, num_points_x)) + return (xpts) +end + +# Cubed sphere +function target_coordinates( + space::Spaces.SpectralElementSpace2D, + num_points, + ::Domains.SphereDomain, +) + num_points_long, num_points_lat = num_points + FT = Spaces.undertype(space) + longpts = collect(range(FT(-180), FT(180), num_points_long)) + latpts = collect(range(FT(-80), FT(80), num_points_lat)) + + return (longpts, latpts) +end + +# Box +function add_space_coordinates_maybe!( + nc::NCDatasets.NCDataset, + space::Spaces.SpectralElementSpace2D, + num_points, + ::Domains.RectangleDomain, +) + xname, yname = ("x", "y") + xpts, ypts = target_coordinates(space, num_points) + add_dimension_maybe!(nc, "x", xpts; units = "m") + add_dimension_maybe!(nc, "y", ypts; units = "m") + return [xname, yname] +end + +# Plane +function add_space_coordinates_maybe!( + nc::NCDatasets.NCDataset, + space::Spaces.SpectralElementSpace1D, + num_points, + ::Domains.IntervalDomain, +) + xname = "x" + xpts = target_coordinates(space, num_points) + add_dimension_maybe!(nc, "x", xpts; units = "m") + return [xname] +end + +# Cubed sphere +function add_space_coordinates_maybe!( + nc::NCDatasets.NCDataset, + space::Spaces.SpectralElementSpace2D, + num_points, + ::Domains.SphereDomain, +) + latname, longname = ("lat", "lon") + latpts, longpts = target_coordinates(space, num_points) + add_dimension_maybe!(nc, "lat", latpts; units = "degrees_north") + add_dimension_maybe!(nc, "lon", longpts; units = "degrees_east") + return [latname, longname] +end + +# General hybrid space. This calls both the vertical and horizontal add_space_coordinates_maybe! +# and combines the resulting dictionaries +function add_space_coordinates_maybe!( + nc::NCDatasets.NCDataset, + space::Spaces.ExtrudedFiniteDifferenceSpace{S}, + num_points, +) where {S <: Spaces.Staggering} + + hdims_names = vdims_names = [] + + num_points_horiz..., num_points_vertic = num_points + + # Being an Extruded space, we can assume that we have an horizontal and a vertical space. + # We can also assume that the vertical space has dimension 1 + hdims_names = add_space_coordinates_maybe!( + nc, + Spaces.horizontal_space(space), + num_points_horiz, + ) + + vertical_space = Spaces.FiniteDifferenceSpace{S}(space.vertical_topology) + vdims_names = + add_space_coordinates_maybe!(nc, vertical_space, num_points_vertic) + + hdims_names == vdims_names == [] && error("Found empty space") + + return vcat(hdims_names, vdims_names) +end + +# General hybrid space. This calls both the vertical and horizontal add_space_coordinates_maybe! +# and combines the resulting dictionaries +function target_coordinates( + space::Spaces.ExtrudedFiniteDifferenceSpace{S}, + num_points, +) where {S <: Spaces.Staggering} + + hcoords = vcoords = () + + num_points_horiz..., num_points_vertic = num_points + + hcoords = + target_coordinates(Spaces.horizontal_space(space), num_points_horiz) + + vertical_space = Spaces.FiniteDifferenceSpace{S}(space.vertical_topology) + vcoords = target_coordinates(vertical_space, num_points_vertic) + + hcoords == vcoords == () && error("Found empty space") + + return hcoords, vcoords +end + +function coordinate_type_from_horizontal_space( + space::Spaces.SpectralElementSpace2D, + domain::Domains.SphereDomain, +) + return Geometry.LatLongPoint +end + + +function coordinate_type_from_horizontal_space( + space::Spaces.SpectralElementSpace2D, + domain::Domains.RectangleDomain, +) + return Geometry.XYPoint +end + +function coordinate_type_from_horizontal_space( + space::Spaces.SpectralElementSpace1D, + domain::Domains.IntervalDomain, +) + return Geometry.XPoint +end + +""" + coordinate_type_from_horizontal_space!(space::Spaces.AbstractSpace) + +Return the `ClimaCore.Geometry` coordinate type for the horizontal space associated to `space`. + +Examples of return values are `ClimaCore.Geometry.LatLongPoint` or `ClimaCore.Geometry.XYPoint`. +""" +coordinate_type_from_horizontal_space(space) = + coordinate_type_from_horizontal_space(space, space.topology.mesh.domain) + +struct NetCDFWriter{T} + + # TODO: At the moment, each variable gets its remapper. This is a little bit of a waste + # because we probably only need a handful of remappers since the same remapper can be + # used for multiple fields as long as they are all defined on the same space. We need + # just a few remappers because realistically we need to support fields defined on the + # entire space and fields defined on 2D slices. However, handling this distinction at + # construction time is quite difficult. + remappers::Dict{String, Remapper} + + # Tuple/Array of integers that identifies how many points to use for interpolation along + # the various dimensions. It has to have the same size as the target interpolation + # space. + num_points::T + + # How much to compress the data in the final NetCDF file: 0 no compression, 9 max + # compression. + compression_level::Int + + # NetCDF files that are currently open. Only the root process uses this field. + open_files::Dict{String, NCDatasets.NCDataset} end +""" + close(writer::NetCDFWriter) + + +Close all the files open in `writer`. +""" +close(writer::NetCDFWriter) = map(close, values(writer.open_files)) """ NetCDFWriter() Save a `ScheduledDiagnostic` to a NetCDF file inside the `output_dir` of the simulation by -performing a pointwise (non-conservative) remapping first. This writer does not work on -distributed simulations. +performing a pointwise (non-conservative) remapping first. +Keyword arguments +================== + +- `num_points`: Lat-Long-Z, X-Y-Z. TODO: This is a very barebone NetCDFWriter. This writer only supports 3D fields on cubed spheres at the moment. Do not consider this implementation as the "final word". @@ -88,70 +444,117 @@ We need to implement the following features/options: - ...more features/options """ -function NetCDFWriter(; - num_points_longitude = 100, - num_points_latitude = 100, - num_points_altitude = 100, - compression_level = 9, -) - function write_to_netcdf(field, diagnostic, integrator) - var = diagnostic.variable - time = integrator.t +function NetCDFWriter(; num_points = (80, 180, 10), compression_level = 9) + return NetCDFWriter{typeof(num_points)}( + Dict(), + num_points, + compression_level, + Dict(), + ) +end - # TODO: Automatically figure out details on the remapping depending on the given - # field +function write_field!(writer::NetCDFWriter, field, diagnostic, integrator) + + var = diagnostic.variable + space = axes(field) + FT = Spaces.undertype(space) + # We have to deal with to cases: when we have an horizontal slice (e.g., the + # surface), and when we have a full space. We distinguish these cases by checking if + # the given space has the horizontal_space attribute. If not, it is going to be a + # SpectralElementSpace2D and we don't have to deal with the z coordinates. + is_horizontal_space = !hasproperty(space, :horizontal_space) + + if !is_horizontal_space + horizontal_space = Spaces.horizontal_space(space) + is_plane = typeof(horizontal_space) <: Spaces.SpectralElementSpace1D + else + horizontal_space = space + is_plane = false + end + # Prepare the remapper if we don't have one for the given variable. We need one remapper + # per variable (not one per diagnostic since all the time reductions return the same + # type of space). + + # TODO: Expand this once we support spatial reductions + if !haskey(writer.remappers, var.short_name) + if is_horizontal_space + hpts = target_coordinates(space, writer.num_points) + vpts = [] + else + hpts, vpts = target_coordinates(space, writer.num_points) + end + + # zcoords is going to be empty for a 2D horizontal slice + zcoords = [Geometry.ZPoint(p) for p in vpts] + + hcoord_type = coordinate_type_from_horizontal_space(horizontal_space) + + if is_plane + hcoords = [hcoord_type(hc1) for hc1 in hpts] + else + hcoords = [hcoord_type(hc1, hc2) for hc1 in hpts[1], hc2 in hpts[2]] + end + + writer.remappers[var.short_name] = Remapper(hcoords, zcoords, space) + end - # Let's support only cubed spheres for the moment - typeof(axes(field).horizontal_space.topology.mesh) <: - Meshes.AbstractCubedSphere || - error("NetCDF writer supports only cubed sphere at the moment") + remapper = writer.remappers[var.short_name] - output_path = joinpath( - integrator.p.simulation.output_dir, - "$(diagnostic.output_short_name)_$time.nc", - ) + # There's an MPI call in here (to aggregate the results) + interpolated_field = interpolate(remapper, field) - vert_domain = axes(field).vertical_topology.mesh.domain - z_min, z_max = vert_domain.coord_min.z, vert_domain.coord_max.z + # Only the root process has to write + ClimaComms.iamroot(ClimaComms.context(field)) || return - FT = Spaces.undertype(axes(field)) + output_path = joinpath( + integrator.p.simulation.output_dir, + "$(diagnostic.output_short_name).nc", + ) - longpts = range( - Geometry.LongPoint(-FT(180.0)), - Geometry.LongPoint(FT(180.0)), - length = num_points_longitude, - ) - latpts = range( - Geometry.LatPoint(-FT(80.0)), - Geometry.LatPoint(FT(80.0)), - length = num_points_latitude, - ) - zpts = range( - Geometry.ZPoint(FT(z_min)), - Geometry.ZPoint(FT(z_max)), - length = num_points_altitude, - ) + if !haskey(writer.open_files, output_path) + # Append or write a new file + open_mode = isfile(output_path) ? "a" : "c" + writer.open_files[output_path] = + NCDatasets.Dataset(output_path, open_mode) + end - nc = NCDatasets.Dataset(output_path, "c") + nc = writer.open_files[output_path] - NCDatasets.defDim(nc, "lon", num_points_longitude) - NCDatasets.defDim(nc, "lat", num_points_latitude) - NCDatasets.defDim(nc, "z", num_points_altitude) + # Define time coordinate + add_time_maybe!(nc, FT) - nc.attrib["long_name"] = diagnostic.output_long_name - nc.attrib["units"] = var.units - nc.attrib["comments"] = var.comments + dim_names = add_space_coordinates_maybe!(nc, space, writer.num_points) + if haskey(nc, "$(var.short_name)") + # We already have something in the file + v = nc["$(var.short_name)"] + spatial_size..., temporal_size = size(v) + spatial_size == size(interpolated_field) || + error("incompatible dimensions for $(var.short_name)") + else v = NCDatasets.defVar( nc, "$(var.short_name)", FT, - ("lon", "lat", "z"), - deflatelevel = compression_level, + (dim_names..., "time"), + deflatelevel = writer.compression_level, ) + v.attrib["long_name"] = diagnostic.output_long_name + v.attrib["units"] = var.units + v.attrib["comments"] = var.comments + + temporal_size = 0 + end + + # We need to write to the next position after what we read from the data (or the first + # position ever if we are writing the file for the first time) + time_index = temporal_size + 1 - v[:, :, :] = interpolate_array(field, longpts, latpts, zpts) + nc["time"][time_index] = integrator.t - close(nc) + if is_horizontal_space || is_plane + v[:, :, time_index] = interpolated_field + else + v[:, :, :, time_index] = interpolated_field end end diff --git a/src/solver/solve.jl b/src/solver/solve.jl index ea08f573eed..9a0e2654b3f 100644 --- a/src/solver/solve.jl +++ b/src/solver/solve.jl @@ -71,6 +71,9 @@ function solve_atmos!(integrator) @error "ClimaAtmos simulation crashed. Stacktrace for failed simulation" exception = (ret_code, catch_backtrace()) return AtmosSolveResults(nothing, :simulation_crashed, nothing) + finally + # Close all the files opened by the writers + map(CAD.close, integrator.p.writers) end end diff --git a/src/solver/type_getters.jl b/src/solver/type_getters.jl index bce90b82aa3..936ba491975 100644 --- a/src/solver/type_getters.jl +++ b/src/solver/type_getters.jl @@ -967,7 +967,12 @@ function get_integrator(config::AtmosConfig) diagnostic_counters[diag] = 1 # If it is not a reduction, call the output writer as well if isnothing(diag.reduction_time_func) - diag.output_writer(diagnostic_storage[diag], diag, integrator) + CAD.write_field!( + diag.output_writer, + diagnostic_storage[diag], + diag, + integrator, + ) else # Add to the accumulator diagnostic_accumulators[diag] =