Pedenatic compiler options and RT instability demo.

CMakeLists.txt

@@ -19,11 +19,11 @@ pybind11_add_module(
 target_link_libraries(${PROJECT_NAME} PUBLIC dolfinx)
 target_compile_options(
   ${PROJECT_NAME} PRIVATE
-  # -Wno-comment
-  # -Wall
-  # -Wextra
-  # -pedantic
-  # -Werror
+  -Wno-comment
+  -Wall
+  -Wextra
+  -pedantic
+  -Werror
   -Wfatal-errors
   )
 

demo/demo_rayleigh-taylor-instability/demo_raleigh-taylor-instability.py

@@ -0,0 +1,165 @@
+# Copyright (c) 2023 Nathan Sime
+# This file is part of LEoPart-X, a particle-in-cell package for DOLFIN-X
+# License: GNU Lesser GPL version 3 or any later version
+# SPDX-License-Identifier:    LGPL-3.0-or-later
+
+from mpi4py import MPI
+from petsc4py import PETSc
+
+import adios2
+import numpy as np
+
+import dolfinx
+import dolfinx.fem.petsc
+import leopart.cpp as pyleopart
+import leopart.io
+import ufl
+
+
+def pprint(*msg, rank=None):
+    if rank is not None and MPI.COMM_WORLD.rank != rank:
+        return
+    print(f"[{MPI.COMM_WORLD.rank}]: {' '.join(map(str, msg))}", flush=True)
+
+
+# Parameters
+lmbda, H = 0.9142, 1.0
+p = 2
+A = 0.02
+db = 0.2
+S = db * (1 - db)
+tableau = pyleopart.tableaus.order2.explicit_midpoint()
+t_max = 2000.0
+num_t_steps_max = 200
+
+# Geometry
+mesh = dolfinx.mesh.create_rectangle(
+    MPI.COMM_WORLD, [[0.0, 0.0], [lmbda, H]], [40, 40],
+    cell_type=dolfinx.mesh.CellType.triangle,
+    diagonal=dolfinx.mesh.DiagonalType.left_right)
+
+# Chemistry space
+PHI = dolfinx.fem.FunctionSpace(mesh, ("DG", 0))
+phi = dolfinx.fem.Function(PHI)
+phi.interpolate(lambda x: np.where(x[1] > db, 1.0, 0.0))
+
+# Perturb the mesh and initial condition
+meshx = mesh.geometry.x
+meshx[:, 1] += meshx[:, 1] * (H - meshx[:, 1]) / S \
+              * A * np.cos(np.pi * lmbda * meshx[:, 0])
+
+# Generate particles and interpolate the composition
+xp, p2cell = pyleopart.mesh_fill(mesh._cpp_object, 30)
+xp = np.c_[xp, np.zeros_like(xp[:, 0])]
+pprint(f"num paticles: {xp.shape[0]}")
+ptcls = pyleopart.Particles(xp, p2cell)
+tableau.check_and_create_fields(ptcls)
+ptcls.add_field("phi", [1])
+pyleopart.transfer_to_particles(ptcls, ptcls.field("phi"), phi._cpp_object)
+
+# Viscosity model
+eta_top = dolfinx.fem.Constant(mesh, 1.0)
+eta_bottom = dolfinx.fem.Constant(mesh, 0.1)
+mu = eta_bottom + phi * (eta_top - eta_bottom)
+
+# Standard Taylor Hood mixed element
+Ve = ufl.VectorElement("CG", mesh.ufl_cell(), p)
+Qe = ufl.FiniteElement("CG", mesh.ufl_cell(), p - 1)
+We = ufl.MixedElement([Ve, Qe])
+W = dolfinx.fem.FunctionSpace(mesh, We)
+u, p_ = ufl.TrialFunctions(W)
+v, q = ufl.TestFunctions(W)
+Uh = dolfinx.fem.Function(W)
+
+# Bilinear formulation
+a = (
+        ufl.inner(2 * mu * ufl.sym(ufl.grad(u)), ufl.grad(v)) * ufl.dx
+        - p_ * ufl.div(v) * ufl.dx
+        - q * ufl.div(u) * ufl.dx
+)
+
+Rb = dolfinx.fem.Constant(mesh, 1.0)
+f = Rb * phi * ufl.as_vector((0, -1))
+L = ufl.inner(f, v) * ufl.dx
+
+# Create BCs: free slip on left and right, zero flow top and bottom
+facets_top_bot = dolfinx.mesh.locate_entities_boundary(
+    mesh, dim=mesh.topology.dim - 1,
+    marker=lambda x: np.isclose(x[1], 0.0) | np.isclose(x[1], H))
+facets_left_right = dolfinx.mesh.locate_entities_boundary(
+    mesh, dim=mesh.topology.dim - 1,
+    marker=lambda x: np.isclose(x[0], 0.0) | np.isclose(x[0], lmbda))
+
+V_x = W.sub(0).sub(0).collapse()
+zero = dolfinx.fem.Function(V_x[0])
+dofs_lr = dolfinx.fem.locate_dofs_topological(
+    (W.sub(0).sub(0), V_x[0]), mesh.topology.dim - 1, facets_left_right)
+zero_x_bc = dolfinx.fem.dirichletbc(zero, dofs_lr, W.sub(0).sub(0))
+
+W0 = W.sub(0).collapse()
+zero = dolfinx.fem.Function(W0[0])
+dofs_tb = dolfinx.fem.locate_dofs_topological(
+    (W.sub(0), W0[0]), mesh.topology.dim - 1, facets_top_bot)
+zero_y_bc = dolfinx.fem.dirichletbc(zero, dofs_tb, W.sub(0))
+
+# Pin pressure DoF in bottom left corner for solvable system
+Q = W.sub(1).collapse()
+zero_p = dolfinx.fem.Function(Q[0])
+dofs_p = dolfinx.fem.locate_dofs_geometrical(
+    (W.sub(1), Q[0]),
+    lambda x: np.isclose(x[0], 0.0) & np.isclose(x[1], 0.0))
+zero_p_bc = dolfinx.fem.dirichletbc(zero_p, dofs_p, W.sub(1))
+
+problem = dolfinx.fem.petsc.LinearProblem(
+    a, L, u=Uh, bcs=[zero_x_bc, zero_y_bc, zero_p_bc],
+    petsc_options={"ksp_type": "preonly", "pc_type": "lu",
+                   "pc_factor_mat_solver_type": "mumps"})
+
+# Velocity as function of time to be used in Runge-Kutta integration
+uh = Uh.sub(0)
+def velocity(t):
+    pyleopart.transfer_to_function(phi._cpp_object, ptcls, ptcls.field("phi"))
+    problem.solve()
+    return uh._cpp_object
+
+# h measured used in CFL criterion estimation
+h_measure = dolfinx.cpp.mesh.h(
+    mesh._cpp_object, mesh.topology.dim - 1, np.arange(
+        mesh.topology.index_map(mesh.topology.dim - 1).size_local,
+        dtype=np.int32))
+hmin = mesh.comm.allreduce(h_measure.min(), op=MPI.MIN)
+
+# Output files
+ptcl_file = leopart.io.XDMFParticlesFile(
+    mesh.comm, "particles.xdmf", adios2.Mode.Write)
+phi_file = dolfinx.io.XDMFFile(mesh.comm, "phi.xdmf", "w")
+
+# Space for estimating speed for CFL criterion
+Vspd = dolfinx.fem.FunctionSpace(mesh, ("DG", 0))
+uspd = dolfinx.fem.Function(Vspd)
+uspd_expr = dolfinx.fem.Expression(
+    ufl.sqrt(ufl.inner(uh, uh)), Vspd.element.interpolation_points())
+
+# Main time loop
+t = 0.0
+ptcl_file.write_particles(ptcls, t, ["phi"])
+phi_file.write_mesh(mesh)
+phi_file.write_function(phi, t=t)
+problem.solve()  # Get initial speed estimate
+for j in range(num_t_steps_max):
+    uspd.interpolate(uspd_expr)
+    uspd.vector.assemble()
+    max_u_vec = uspd.vector.norm(PETSc.NormType.INF)
+    dt = (c_cfl := 1.0) * hmin / max_u_vec
+    t += dt
+    pprint(f"Time step {j}, dt = {dt:.3e}, t = {t:.3e}", rank=0)
+
+    pyleopart.rk(mesh._cpp_object, ptcls, tableau,
+                 velocity, t, dt)
+
+    pyleopart.transfer_to_function(phi._cpp_object, ptcls, ptcls.field("phi"))
+    ptcl_file.write_particles(ptcls, t, ["phi"])
+    phi_file.write_function(phi, t=t)
+
+    if t > t_max:
+        break

leopart/cpp/Particles.h

@@ -9,6 +9,7 @@
 #include <vector>
 #include <stdexcept>
 #include <map>
+#include <ranges>
 
 
 namespace dolfinx
@@ -94,16 +95,15 @@ class Particles
   /// data.
   inline std::vector<std::size_t> active_pidxs()
   {
-    const std::size_t num_pidxs = _particle_to_cell.size();
-    std::vector<std::size_t> valid_pidxs(num_pidxs);
-    std::size_t n_valid_particles = 0;
-    for (std::size_t pidx = 0; pidx < num_pidxs; ++pidx)
-    {
-      if (_particle_to_cell[pidx] != INVALID_CELL)
-        valid_pidxs[n_valid_particles++] = pidx;
-    }
-    valid_pidxs.resize(n_valid_particles);
-    return valid_pidxs;
+    const auto is_pidx_valid =
+      [&invalid_cell=INVALID_CELL, &p2c=_particle_to_cell](
+        const std::size_t& p_idx)
+      {
+        return p2c[p_idx] != invalid_cell;
+      };
+    auto v = std::views::iota((std::size_t) 0, _particle_to_cell.size())
+           | std::views::filter(is_pidx_valid);
+    return std::vector<std::size_t>(v.begin(), v.end());
   }
 
   /// Given a process local particle index, return the owning cell and cell

leopart/cpp/advect.h

@@ -260,7 +260,7 @@ void rk(
   dolfinx::common::Timer timer("leopart::advect::rk");
 
   const std::string xn_name = tableau.field_name_xn();
-  const int num_steps = tableau.order;
+  const std::size_t num_steps = tableau.order;
   mdspan_t<const T, 2> a = tableau.a_md();
   const std::vector<T>& b = tableau.b;
   const std::vector<T>& c = tableau.c;

leopart/cpp/generation.h

@@ -63,7 +63,6 @@ mesh_fill(
 
   // Geometric dimension
   const std::size_t gdim = mesh.geometry().dim();
-  const int tdim = mesh.topology()->dim();
 
   // Get coordinate map
   if (mesh.geometry().cmaps().size() > 1)
@@ -195,7 +194,7 @@ std::vector<T> random_reference_tetrahedron(
   std::mt19937 rgen(seed);
   std::uniform_real_distribution<T> dist(-1.0, 1.0);
 
-  for (int i = 0; i < n; ++i)
+  for (std::size_t i = 0; i < n; ++i)
   {
     std::generate(r.begin(), r.end(),
                   [&dist, &rgen]() { return dist(rgen); });
@@ -222,7 +221,7 @@ std::vector<T> random_reference_tetrahedron(
     }
 
     const auto p_row = &p[i * gdim];
-    for (int j = 0; j < gdim; ++j)
+    for (std::size_t j = 0; j < gdim; ++j)
       p_row[j] += r[j];
   }
 

leopart/cpp/transfer.h

@@ -88,8 +88,8 @@ void transfer_to_function_l2_callback(
   // @todo these definitions are legacy DOLFIN and should be refactored for
   // appropriate unrolling of DoFs.
   const int block_size = element->block_size();
-  const int value_size = element->value_size() / block_size;
-  const int space_dimension = element->space_dimension() / block_size;
+  // const int value_size = element->value_size() / block_size;
+  const std::size_t space_dimension = element->space_dimension() / block_size;
 
   std::shared_ptr<const dolfinx::fem::DofMap> dm
       = f->function_space()->dofmap();
@@ -115,10 +115,10 @@ void transfer_to_function_l2_callback(
   //  QT . Q: (space_dimension x space_dimension)
   //  L: (n_p x block_size)
   //  QT . L: (space_dimension x block_size)
-  for (int c = 0; c < ncells; ++c)
+  for (std::int32_t c = 0; c < ncells; ++c)
   {
     const std::vector<std::size_t> cell_particles = cell_to_particle[c];
-    int cell_np = cell_particles.size();
+    const std::size_t cell_np = cell_particles.size();
 
     // Assemble Q
     std::vector<T> Q_T_data(cell_np * space_dimension);
@@ -129,7 +129,7 @@ void transfer_to_function_l2_callback(
     for (std::size_t cell_p = 0; cell_p < cell_np; ++cell_p)
     {
       const std::size_t p_idx = cell_particles[cell_p];
-      for (int i = 0; i < space_dimension; ++i)
+      for (std::size_t i = 0; i < space_dimension; ++i)
         Q(cell_p, i) = basis_evals_md(p_idx, i, 0); // Assume no vector valued basis
     }
     leopart::math::transpose<T>(Q, Q_T);
@@ -157,7 +157,7 @@ void transfer_to_function_l2_callback(
 
     // Populate FE function DoFs
     const auto& dofs = dm->cell_dofs(c);
-    for (int i = 0; i < dofs.size(); ++i)
+    for (std::size_t i = 0; i < dofs.size(); ++i)
       for (int k = 0; k < block_size; ++k)
         expansion_coefficients[dofs[i]*block_size + k] = soln[i*block_size + k];
   }
@@ -230,7 +230,7 @@ void transfer_to_function_constrained(
                              space_dimension * value_size * 2);
   std::vector<T> ce0(0, 0.0), ci0(space_dimension * value_size * 2, 0.0);
 
-  for (std::size_t i = 0; i < space_dimension; i++)
+  for (int i = 0; i < space_dimension; i++)
   {
     CI(i, i) = 1.;
     CI(i, i + space_dimension) = -1;

leopart/cpp/wrapper.cpp

@@ -61,7 +61,7 @@ PYBIND11_MODULE(cpp, m)
       .def(py::init(
             [](const py::array_t<dtype, py::array::c_style>& px,
                const std::vector<std::int32_t>& p_cells) {
-                if (px.shape(0) != p_cells.size())
+                if ((std::size_t) px.shape(0) != p_cells.size())
                   throw std::invalid_argument(
                     "Number of particles and particle cells should be equivalent");
 
@@ -105,6 +105,13 @@ PYBIND11_MODULE(cpp, m)
               const std::size_t np_per_cell) {
             self.generate_minimum_particles_per_cell(mesh, np_per_cell);
            })
+      .def("active_pidxs",
+           [](Particles<dtype>& self) {
+            const std::vector<std::size_t> pidxs = self.active_pidxs();
+            return py::array_t<std::size_t, py::array::c_style>(
+              pidxs.size(), pidxs.data());
+           },
+           py::return_value_policy::move)
       .def("field_exists", &Particles<dtype>::Particles::field_exists);
 
   // Generation functions

leopart/io.py

@@ -99,11 +99,12 @@ def write_particles(
             field_names: typing.Optional[typing.Sequence[str]] = None):
         if field_names is None:
             field_names = []
+        active_pidxs = particles.active_pidxs()
         data_map = dict(
             zip(field_names,
-                (particles.field(field_name).data()
+                (particles.field(field_name).data()[active_pidxs]
                  for field_name in field_names)))
-        self.write_points(particles.x().data(), data_map, t)
+        self.write_points(particles.x().data()[active_pidxs], data_map, t)
 
     def write_points(
             self, points: np.typing.NDArray,