diff --git a/cpp/demo/interpolation-io/main.cpp b/cpp/demo/interpolation-io/main.cpp index 54a45180d27..7342de4a483 100644 --- a/cpp/demo/interpolation-io/main.cpp +++ b/cpp/demo/interpolation-io/main.cpp @@ -65,7 +65,7 @@ void interpolate_scalar(std::shared_ptr> mesh, #endif } -// This function interpolations a function is a H(curl) finite element +// This function interpolates a function in a H(curl) finite element // space. To visualise the function, it interpolates the H(curl) finite // element function in a discontinuous Lagrange space and outputs the // Lagrange finite element function to a VTX file for visualisation. diff --git a/cpp/dolfinx/fem/utils.h b/cpp/dolfinx/fem/utils.h index db6d5a4332c..dd89a92242f 100644 --- a/cpp/dolfinx/fem/utils.h +++ b/cpp/dolfinx/fem/utils.h @@ -709,8 +709,8 @@ FunctionSpace create_functionspace( /// @brief Create a FunctionSpace from UFC data. /// @param[in] fptr Pointer to a ufcx_function_space_create function. -/// @param[in] function_name Name of a function whose function space to -/// create. Function name is the name of Python variable for +/// @param[in] function_name Name of a function whose function space is to +/// create. Function name is the name of the Python variable for /// ufl.Coefficient, ufl.TrialFunction or ufl.TestFunction as defined in /// the UFL file. /// @param[in] mesh Mesh diff --git a/cpp/dolfinx/io/ADIOS2Writers.h b/cpp/dolfinx/io/ADIOS2Writers.h index 3c0a7464a4f..58169987643 100644 --- a/cpp/dolfinx/io/ADIOS2Writers.h +++ b/cpp/dolfinx/io/ADIOS2Writers.h @@ -780,8 +780,9 @@ void vtx_write_mesh(adios2::IO& io, adios2::Engine& engine, /// @param[in] engine The ADIOS2 engine object /// @param[in] V The function space template -void vtx_write_mesh_from_space(adios2::IO& io, adios2::Engine& engine, - const fem::FunctionSpace& V) +std::pair, std::vector> +vtx_write_mesh_from_space(adios2::IO& io, adios2::Engine& engine, + const fem::FunctionSpace& V) { auto mesh = V.mesh(); assert(mesh); @@ -839,9 +840,18 @@ void vtx_write_mesh_from_space(adios2::IO& io, adios2::Engine& engine, engine.Put(ghost, x_ghost.data()); engine.PerformPuts(); + return {std::move(x_id), std::move(x_ghost)}; } } // namespace impl_vtx +/// Mesh reuse policy +enum class VTXMeshPolicy +{ + update, ///< Re-write the mesh to file upon every write of a fem::Function + reuse ///< Write the mesh to file only the first time a fem::Function is + ///< written to file +}; + /// @brief Writer for meshes and functions using the ADIOS2 VTX format, /// see /// https://adios2.readthedocs.io/en/latest/ecosystem/visualization.html#using-vtk-and-paraview. @@ -865,7 +875,8 @@ class VTXWriter : public ADIOS2Writer VTXWriter(MPI_Comm comm, const std::filesystem::path& filename, std::shared_ptr> mesh, std::string engine = "BPFile") - : ADIOS2Writer(comm, filename, "VTX mesh writer", engine), _mesh(mesh) + : ADIOS2Writer(comm, filename, "VTX mesh writer", engine), _mesh(mesh), + _mesh_reuse_policy(VTXMeshPolicy::update) { // Define VTK scheme attribute for mesh std::string vtk_scheme = impl_vtx::create_vtk_schema({}, {}).str(); @@ -883,12 +894,16 @@ class VTXWriter : public ADIOS2Writer /// element family and degree, and degree-of-freedom map (up to the /// blocksize) must be the same for all functions. /// @param[in] engine ADIOS2 engine type. + /// @param[in] mesh_policy Controls if the mesh is written to file at + /// the first time step only or is re-written (updated) at each time + /// step. /// @note This format supports arbitrary degree meshes. VTXWriter(MPI_Comm comm, const std::filesystem::path& filename, - const typename adios2_writer::U& u, - std::string engine = "BPFile") + const typename adios2_writer::U& u, std::string engine, + VTXMeshPolicy mesh_policy = VTXMeshPolicy::update) : ADIOS2Writer(comm, filename, "VTX function writer", engine), - _mesh(impl_adios2::extract_common_mesh(u)), _u(u) + _mesh(impl_adios2::extract_common_mesh(u)), + _mesh_reuse_policy(mesh_policy), _u(u) { if (u.empty()) throw std::runtime_error("VTXWriter fem::Function list is empty."); @@ -960,6 +975,27 @@ class VTXWriter : public ADIOS2Writer impl_adios2::define_attribute(*_io, "vtk.xml", vtk_scheme); } + /// @brief Create a VTX writer for a list of fem::Functions using + /// the default ADIOS2 engine. + /// + /// This format supports arbitrary degree meshes. + /// + /// @param[in] comm The MPI communicator to open the file on + /// @param[in] filename Name of output file + /// @param[in] u List of functions. The functions must (1) share the + /// same mesh and (2) be (discontinuous) Lagrange functions. The + /// element family and degree must be the same for all functions. + /// @param[in] mesh_policy Controls if the mesh is written to file at + /// the first time step only or is re-written (updated) at each time + /// step. + /// @note This format supports arbitrary degree meshes. + VTXWriter(MPI_Comm comm, const std::filesystem::path& filename, + const typename adios2_writer::U& u, + VTXMeshPolicy mesh_policy = VTXMeshPolicy::update) + : VTXWriter(comm, filename, u, "BPFile", mesh_policy) + { + } + // Copy constructor VTXWriter(const VTXWriter&) = delete; @@ -975,14 +1011,15 @@ class VTXWriter : public ADIOS2Writer // Copy assignment VTXWriter& operator=(const VTXWriter&) = delete; - /// @brief Write data with a given time - /// @param[in] t The time step + /// @brief Write data with a given time stamp. + /// @param[in] t Time stamp to associate with output. void write(double t) { assert(_io); - assert(_engine); adios2::Variable var_step = impl_adios2::define_variable(*_io, "step"); + + assert(_engine); _engine->BeginStep(); _engine->template Put(var_step, t); @@ -991,18 +1028,37 @@ class VTXWriter : public ADIOS2Writer impl_vtx::vtx_write_mesh(*_io, *_engine, *_mesh); else { - // Write a single mesh for functions as they share finite element - std::visit( - [&](auto& u) { - impl_vtx::vtx_write_mesh_from_space(*_io, *_engine, - *u->function_space()); - }, - _u[0]); + if (_mesh_reuse_policy == VTXMeshPolicy::update + or !(_io->template InquireVariable("connectivity"))) + { + // Write a single mesh for functions as they share finite + // element + std::tie(_x_id, _x_ghost) = std::visit( + [&](auto& u) + { + return impl_vtx::vtx_write_mesh_from_space(*_io, *_engine, + *u->function_space()); + }, + _u[0]); + } + else + { + // Node global ids + adios2::Variable orig_id = impl_adios2::define_variable( + *_io, "vtkOriginalPointIds", {}, {}, {_x_id.size()}); + _engine->Put(orig_id, _x_id.data()); + adios2::Variable ghost = impl_adios2::define_variable( + *_io, "vtkGhostType", {}, {}, {_x_ghost.size()}); + _engine->Put(ghost, _x_ghost.data()); + _engine->PerformPuts(); + } // Write function data for each function to file for (auto& v : _u) + { std::visit( [&](auto& u) { impl_vtx::vtx_write_data(*_io, *_engine, *u); }, v); + } } _engine->EndStep(); @@ -1011,6 +1067,12 @@ class VTXWriter : public ADIOS2Writer private: std::shared_ptr> _mesh; adios2_writer::U _u; + + // Control whether the mesh is written to file once or at every time + // step + VTXMeshPolicy _mesh_reuse_policy; + std::vector _x_id; + std::vector _x_ghost; }; /// Type deduction diff --git a/cpp/test/io.cpp b/cpp/test/io.cpp index 3175ee1384e..985794fb6f5 100644 --- a/cpp/test/io.cpp +++ b/cpp/test/io.cpp @@ -77,6 +77,39 @@ void test_fides_function() writer1.write(0.0); } +template +void test_vtx_reuse_mesh() +{ + auto mesh = std::make_shared>( + mesh::create_rectangle(MPI_COMM_WORLD, {{{0.0, 0.0}, {1.0, 1.0}}}, + {22, 12}, mesh::CellType::triangle)); + + // Create a Basix continuous Lagrange element of degree 1 + basix::FiniteElement e = basix::create_element( + basix::element::family::P, + mesh::cell_type_to_basix_type(mesh::CellType::triangle), 1, + basix::element::lagrange_variant::unset, + basix::element::dpc_variant::unset, false); + + // Create a scalar function space + auto V = std::make_shared>( + fem::create_functionspace(mesh, e)); + + // Create a finite element Function + auto u = std::make_shared>(V); + auto v = std::make_shared>>(V); + + std::filesystem::path f = "test_vtx_reuse_mesh" + std::to_string(sizeof(T)) + + ".bp"; + io::VTXWriter writer(mesh->comm(), f, {u, v}, "BPFile", + io::VTXMeshPolicy::reuse); + writer.write(0); + + std::fill(u->x()->mutable_array().begin(), u->x()->mutable_array().end(), 1); + + writer.write(1); + +} } // namespace TEST_CASE("Fides output") @@ -87,4 +120,10 @@ TEST_CASE("Fides output") CHECK_NOTHROW(test_fides_function()); } +TEST_CASE("VTX reuse mesh") +{ + CHECK_NOTHROW(test_vtx_reuse_mesh()); + CHECK_NOTHROW(test_vtx_reuse_mesh()); +} + #endif diff --git a/python/dolfinx/fem/petsc.py b/python/dolfinx/fem/petsc.py index 96fdf54defb..078a8338fee 100644 --- a/python/dolfinx/fem/petsc.py +++ b/python/dolfinx/fem/petsc.py @@ -379,7 +379,7 @@ def assemble_matrix_mat(A: PETSc.Mat, a: Form, bcs: list[DirichletBC] = [], def assemble_matrix_nest(a: list[list[Form]], bcs: list[DirichletBC] = [], mat_types=[], diagonal: float = 1.0, constants=None, coeffs=None) -> PETSc.Mat: - """Create a nested matrix and assembled bilinear forms into the matrix. + """Create a nested matrix and assemble bilinear forms into the matrix. Args: a: Rectangular (list-of-lists) array for bilinear forms. diff --git a/python/dolfinx/io/__init__.py b/python/dolfinx/io/__init__.py index 88d90108d6a..5804ddd3978 100644 --- a/python/dolfinx/io/__init__.py +++ b/python/dolfinx/io/__init__.py @@ -13,5 +13,6 @@ if _cpp.common.has_adios2: # FidesWriter and VTXWriter require ADIOS2 - from dolfinx.io.utils import FidesMeshPolicy, FidesWriter, VTXWriter - __all__ = [*__all__, "FidesWriter", "VTXWriter", "FidesMeshPolicy"] + from dolfinx.io.utils import FidesMeshPolicy, FidesWriter, VTXMeshPolicy, VTXWriter + __all__ = [*__all__, "FidesWriter", "VTXWriter", "FidesMeshPolicy", + "VTXMeshPolicy"] diff --git a/python/dolfinx/io/utils.py b/python/dolfinx/io/utils.py index 5d342443f29..0fca232ec79 100644 --- a/python/dolfinx/io/utils.py +++ b/python/dolfinx/io/utils.py @@ -37,8 +37,9 @@ def _extract_cpp_functions(functions: typing.Union[list[Function], Function]): # FidesWriter and VTXWriter require ADIOS2 if _cpp.common.has_adios2: - from dolfinx.cpp.io import FidesMeshPolicy # F401 - __all__ = [*__all__, "FidesWriter", "VTXWriter", "FidesMeshPolicy"] + from dolfinx.cpp.io import FidesMeshPolicy, VTXMeshPolicy # F401 + __all__ = [*__all__, "FidesWriter", "VTXWriter", "FidesMeshPolicy", + "VTXMeshPolicy"] class VTXWriter: """Writer for VTX files, using ADIOS2 to create the files. @@ -54,7 +55,8 @@ class VTXWriter: def __init__(self, comm: _MPI.Comm, filename: typing.Union[str, Path], output: typing.Union[Mesh, Function, list[Function]], - engine: str = "BPFile"): + engine: str = "BPFile", + mesh_policy: VTXMeshPolicy = VTXMeshPolicy.update): """Initialize a writer for outputting data in the VTX format. Args: @@ -65,6 +67,11 @@ def __init__(self, comm: _MPI.Comm, filename: typing.Union[str, Path], (discontinuous Lagrange Functions. engine: ADIOS2 engine to use for output. See ADIOS2 documentation for options. + mesh_policy: Controls if the mesh is written to file at + the first time step only when a ``Function`` is + written to file, or is re-written (updated) at each + time step. Has an effect only for ``Function`` + output. Note: All Functions for output must share the same mesh and @@ -90,7 +97,7 @@ def __init__(self, comm: _MPI.Comm, filename: typing.Union[str, Path], except (NotImplementedError, TypeError, AttributeError): # Input is a single function or a list of functions self._cpp_object = _vtxwriter(comm, filename, _extract_cpp_functions( - output), engine) # type: ignore[arg-type] + output), engine, mesh_policy) # type: ignore[arg-type] def __enter__(self): return self diff --git a/python/dolfinx/la.py b/python/dolfinx/la.py index 9d2702233a6..97cc7ae8d63 100644 --- a/python/dolfinx/la.py +++ b/python/dolfinx/la.py @@ -214,7 +214,7 @@ def norm(self, type=_cpp.la.Norm.l2) -> np.floating: """Compute a norm of the vector. Args: - type: Norm type to computed. + type: Norm type to compute. Returns: Computed norm. diff --git a/python/dolfinx/wrappers/io.cpp b/python/dolfinx/wrappers/io.cpp index 3c0bbf3e6ab..ac07f9f66bb 100644 --- a/python/dolfinx/wrappers/io.cpp +++ b/python/dolfinx/wrappers/io.cpp @@ -130,12 +130,13 @@ void declare_vtx_writer(nb::module_& m, std::string type) const dolfinx::fem::Function, T>>, std::shared_ptr, T>>>>& u, - std::string engine) { + std::string engine, dolfinx::io::VTXMeshPolicy policy) { new (self) - dolfinx::io::VTXWriter(comm.get(), filename, u, engine); + dolfinx::io::VTXWriter(comm.get(), filename, u, engine, policy); }, nb::arg("comm"), nb::arg("filename"), nb::arg("u"), - nb::arg("engine")) + nb::arg("engine") = "BPFile", + nb::arg("policy") = dolfinx::io::VTXMeshPolicy::update) .def("close", [](dolfinx::io::VTXWriter& self) { self.close(); }) .def( "write", @@ -338,6 +339,10 @@ void io(nb::module_& m) nb::enum_(m, "FidesMeshPolicy") .value("update", dolfinx::io::FidesMeshPolicy::update) .value("reuse", dolfinx::io::FidesMeshPolicy::reuse); + + nb::enum_(m, "VTXMeshPolicy") + .value("update", dolfinx::io::VTXMeshPolicy::update) + .value("reuse", dolfinx::io::VTXMeshPolicy::reuse); #endif declare_vtx_writer(m, "float32"); diff --git a/python/test/unit/io/test_adios2.py b/python/test/unit/io/test_adios2.py index 22cee7c62ac..a48c0b203c3 100644 --- a/python/test/unit/io/test_adios2.py +++ b/python/test/unit/io/test_adios2.py @@ -20,9 +20,9 @@ from dolfinx.mesh import CellType, create_mesh, create_unit_cube, create_unit_square try: - from dolfinx.io import FidesWriter, VTXWriter + from dolfinx.io import FidesWriter, VTXMeshPolicy, VTXWriter except ImportError: - pytest.skip("Test require ADIOS2", allow_module_level=True) + pytest.skip("Tests require ADIOS2", allow_module_level=True) def generate_mesh(dim: int, simplex: bool, N: int = 5, dtype=None): @@ -84,7 +84,7 @@ def vel(x): @pytest.mark.skipif(not has_adios2, reason="Requires ADIOS2.") @pytest.mark.parametrize("dim", [2, 3]) @pytest.mark.parametrize("simplex", [True, False]) -def test_findes_single_function(tempdir, dim, simplex): +def test_fides_single_function(tempdir, dim, simplex): "Test saving a single first order Lagrange functions" mesh = generate_mesh(dim, simplex) v = Function(functionspace(mesh, ("Lagrange", 1))) @@ -288,3 +288,39 @@ def partitioner(comm, nparts, local_graph, num_ghost_nodes): filename = Path(tempdir, "empty_rank_mesh.bp") with VTXWriter(comm, filename, u) as f: f.write(0.0) + + +@pytest.mark.skipif(not has_adios2, reason="Requires ADIOS2.") +@pytest.mark.parametrize("dim", [2, 3]) +@pytest.mark.parametrize("simplex", [True, False]) +@pytest.mark.parametrize("reuse", [True, False]) +def test_vtx_reuse_mesh(tempdir, dim, simplex, reuse): + "Test reusage of mesh by VTXWriter" + adios2 = pytest.importorskip("adios2") + + mesh = generate_mesh(dim, simplex) + v = Function(functionspace(mesh, ("Lagrange", 1))) + filename = Path(tempdir, "v.bp") + v.name = "v" + policy = VTXMeshPolicy.reuse if reuse else VTXMeshPolicy.update + + # Save three steps + writer = VTXWriter(mesh.comm, filename, v, "BP4", policy) + writer.write(0) + v.interpolate(lambda x: 0.5 * x[0]) + writer.write(1) + v.interpolate(lambda x: x[1]) + writer.write(2) + writer.close() + + reuse_variables = ["NumberOfEntities", "NumberOfNodes", "connectivity", "geometry", "types"] + target_all = 3 # For all other variables the step count is number of writes + target_mesh = 1 if reuse else 3 # For mesh variables the step count is 1 if reuse else number of writes + + adios_file = adios2.open(str(filename), "r", comm=mesh.comm, engine_type="BP4") + for name, var in adios_file.available_variables().items(): + if name in reuse_variables: + assert int(var["AvailableStepsCount"]) == target_mesh + else: + assert int(var["AvailableStepsCount"]) == target_all + adios_file.close()