From bc6097c76bdf7397d51a1a76a0ed25d76e73f90c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= Date: Wed, 21 Aug 2024 11:53:13 +0000 Subject: [PATCH 1/4] Skip negative cell indices when packing coefficients --- cpp/dolfinx/fem/utils.h | 9 + python/test/unit/fem/test_assemble_submesh.py | 169 +++++++++++++++++- 2 files changed, 177 insertions(+), 1 deletion(-) diff --git a/cpp/dolfinx/fem/utils.h b/cpp/dolfinx/fem/utils.h index 619a9e02d7c..77287e3953b 100644 --- a/cpp/dolfinx/fem/utils.h +++ b/cpp/dolfinx/fem/utils.h @@ -914,6 +914,8 @@ void pack_coefficient_entity(std::span c, int cstride, { auto entity = entities.subspan(e, estride); std::int32_t cell = fetch_cells(entity); + if (cell < 0) + continue; auto cell_coeff = c.subspan((e / estride) * cstride + offset, space_dim); pack(cell_coeff, cell, bs, v, cell_info, dofmap, transformation); } @@ -923,6 +925,8 @@ void pack_coefficient_entity(std::span c, int cstride, { auto entity = entities.subspan(e, estride); std::int32_t cell = fetch_cells(entity); + if (cell < 0) + continue; auto cell_coeff = c.subspan((e / estride) * cstride + offset, space_dim); pack(cell_coeff, cell, bs, v, cell_info, dofmap, transformation); } @@ -932,6 +936,8 @@ void pack_coefficient_entity(std::span c, int cstride, { auto entity = entities.subspan(e, estride); std::int32_t cell = fetch_cells(entity); + if (cell < 0) + continue; auto cell_coeff = c.subspan(e / estride * cstride + offset, space_dim); pack(cell_coeff, cell, bs, v, cell_info, dofmap, transformation); } @@ -941,6 +947,8 @@ void pack_coefficient_entity(std::span c, int cstride, { auto entity = entities.subspan(e, estride); std::int32_t cell = fetch_cells(entity); + if (cell < 0) + continue; auto cell_coeff = c.subspan((e / estride) * cstride + offset, space_dim); pack(cell_coeff, cell, bs, v, cell_info, dofmap, transformation); } @@ -1115,6 +1123,7 @@ void pack_coefficients(const Form& form, IntegralType integral_type, auto mesh = coefficients[coeff]->function_space()->mesh(); std::vector facets = form.domain(IntegralType::interior_facet, id, *mesh); + std::span cell_info = impl::get_cell_orientation_info(*coefficients[coeff]); diff --git a/python/test/unit/fem/test_assemble_submesh.py b/python/test/unit/fem/test_assemble_submesh.py index 7f69a433f21..15e68bb042d 100644 --- a/python/test/unit/fem/test_assemble_submesh.py +++ b/python/test/unit/fem/test_assemble_submesh.py @@ -1,4 +1,4 @@ -# Copyright (C) 2022 Joseph P. Dean +# Copyright (C) 2022-2024 Joseph P. Dean, Jørgen S. Dokken # # This file is part of DOLFINx (https://www.fenicsproject.org) # @@ -13,13 +13,17 @@ import ufl from dolfinx import default_scalar_type, fem, la +from dolfinx.cpp.fem import compute_integration_domains from dolfinx.mesh import ( GhostMode, + compute_incident_entities, create_box, create_rectangle, create_submesh, create_unit_cube, + create_unit_interval, create_unit_square, + entities_to_geometry, exterior_facet_indices, locate_entities, locate_entities_boundary, @@ -384,3 +388,166 @@ def coeff_expr(x): # TODO Test random mesh and interior facets + + +def test_disjoint_submeshes(): + """Test assembly with multiple disjoint submeshes in same variational form""" + + N = 10 + tol = 1e-14 + mesh = create_unit_interval(MPI.COMM_WORLD, N, ghost_mode=GhostMode.shared_facet) + tdim = mesh.topology.dim + dx = 1.0 / N + center_tag = 1 + left_tag = 2 + right_tag = 3 + left_interface_tag = 4 + right_interface_tag = 5 + + def left(x): + return x[0] < N // 3 * dx + tol + + def right(x): + return x[0] > 2 * N // 3 * dx - tol + + cell_map = mesh.topology.index_map(tdim) + num_cells_local = cell_map.size_local + cell_map.num_ghosts + values = np.full(num_cells_local, center_tag, dtype=np.int32) + values[locate_entities(mesh, tdim, left)] = left_tag + values[locate_entities(mesh, tdim, right)] = right_tag + + cell_tag = meshtags(mesh, tdim, np.arange(num_cells_local, dtype=np.int32), values) + left_facets = compute_incident_entities(mesh.topology, cell_tag.find(left_tag), tdim, tdim - 1) + center_facets = compute_incident_entities( + mesh.topology, cell_tag.find(center_tag), tdim, tdim - 1 + ) + right_facets = compute_incident_entities( + mesh.topology, cell_tag.find(right_tag), tdim, tdim - 1 + ) + + # Create parent facet tag where left interface is tagged with 4, right with 5 + left_interface = np.intersect1d(left_facets, center_facets) + right_interface = np.intersect1d(right_facets, center_facets) + facet_map = mesh.topology.index_map(tdim) + num_facet_local = facet_map.size_local + cell_map.num_ghosts + facet_values = np.full(num_facet_local, 1, dtype=np.int32) + facet_values[left_interface] = left_interface_tag + facet_values[right_interface] = right_interface_tag + facet_tag = meshtags(mesh, tdim - 1, np.arange(num_facet_local, dtype=np.int32), facet_values) + + # Create facet integrals on each interface + left_mesh, left_to_parent, _, _ = create_submesh(mesh, tdim, cell_tag.find(left_tag)) + right_mesh, right_to_parent, _, _ = create_submesh(mesh, tdim, cell_tag.find(right_tag)) + + # One sided interface integral uses only "+" restriction. Sort integration entities such that + # this is always satisfied + def compute_mapped_interior_facet_data(mesh, facet_tag, value, parent_to_sub_map): + """ + Compute integration data for interior facet integrals, where the positive restriction is + always taken on the side that has a cell in the sub mesh + + Parameters + mesh: Parent mesh + facet_tag: Meshtags object for facets + value: Value of the facets to extract + parent_to_sub_map: Mapping from parent mesh to sub mesh + Returns + integration_data: Integration data for interior facets + """ + mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim) + integration_data = compute_integration_domains( + fem.IntegralType.interior_facet, mesh.topology, facet_tag.find(value), facet_tag.dim + ) + mapped_cell_0 = parent_to_sub_map[integration_data[0::4]] + mapped_cell_1 = parent_to_sub_map[integration_data[2::4]] + switch = mapped_cell_1 > mapped_cell_0 + # Order restriction on one side + ordered_integration_data = integration_data.reshape(-1, 4).copy() + if True in switch: + ordered_integration_data[switch, [0, 1, 2, 3]] = ordered_integration_data[ + switch, [2, 3, 0, 1] + ] + return (value, ordered_integration_data.reshape(-1)) + + parent_to_left = np.full(num_cells_local, -1, dtype=np.int32) + parent_to_right = np.full(num_cells_local, -1, dtype=np.int32) + parent_to_left[left_to_parent] = np.arange(len(left_to_parent)) + parent_to_right[right_to_parent] = np.arange(len(right_to_parent)) + integral_data = [ + compute_mapped_interior_facet_data(mesh, facet_tag, left_interface_tag, parent_to_left), + compute_mapped_interior_facet_data(mesh, facet_tag, right_interface_tag, parent_to_right), + ] + + dS = ufl.Measure("dS", domain=mesh, subdomain_data=integral_data) + + def f_left(x): + return np.sin(x[0]) + + def f_right(x): + return x[0] + + V_left = fem.functionspace(left_mesh, ("Lagrange", 1)) + u_left = fem.Function(V_left) + u_left.interpolate(f_left) + + V_right = fem.functionspace(right_mesh, ("Lagrange", 1)) + u_right = fem.Function(V_right) + u_right.interpolate(f_right) + + # Create single integral with different submeshes restrictions + x = ufl.SpatialCoordinate(mesh) + res = "+" + J = x[0] * u_left(res) * dS(left_interface_tag) + ufl.cos(x[0]) * u_right(res) * dS( + right_interface_tag + ) + + # We create an entity map from the parent mesh to the submesh, where + # the cell on either side of the interface is mapped to the same cell + mesh.topology.create_connectivity(tdim - 1, tdim) + f_to_c = mesh.topology.connectivity(tdim - 1, tdim) + parent_to_left = np.full(num_cells_local, -1, dtype=np.int32) + parent_to_right = np.full(num_cells_local, -1, dtype=np.int32) + parent_to_left[left_to_parent] = np.arange(len(left_to_parent)) + parent_to_right[right_to_parent] = np.arange(len(right_to_parent)) + for tag in [4, 5]: + for facet in facet_tag.find(tag): + cells = f_to_c.links(facet) + assert len(cells) == 2 + left_map = parent_to_left[cells] + right_map = parent_to_right[cells] + parent_to_left[cells] = max(left_map) + parent_to_right[cells] = max(right_map) + + entity_maps = {left_mesh: parent_to_left, right_mesh: parent_to_right} + J_compiled = fem.form(J, entity_maps=entity_maps) + J_local = fem.assemble_scalar(J_compiled) + J_sum = mesh.comm.allreduce(J_local, op=MPI.SUM) + + # Compute value of expression at left interface + if len(facets := facet_tag.find(left_interface_tag)) > 0: + left_vertex = entities_to_geometry(mesh, mesh.topology.dim - 1, facets) + left_coord = mesh.geometry.x[left_vertex].reshape(3, -1) + left_val = left_coord[0, 0] * f_left(left_coord)[0] + else: + left_val = 0 + + # Compute value of expression at right interface + if len(facets := facet_tag.find(right_interface_tag)) > 0: + right_vertex = entities_to_geometry(mesh, mesh.topology.dim - 1, facets) + right_coord = mesh.geometry.x[right_vertex].reshape(3, -1) + right_val = np.cos(right_coord[0, 0]) * f_right(right_coord)[0] + else: + right_val = 0 + + # As multiple process might have the vertex, only get value from one process + _, left_rank = mesh.comm.allreduce([abs(left_val), mesh.comm.rank], op=MPI.MAXLOC) + _, right_rank = mesh.comm.allreduce([abs(right_val), mesh.comm.rank], op=MPI.MAXLOC) + if mesh.comm.rank == left_rank: + mesh.comm.send(left_val, dest=0, tag=0) + if mesh.comm.rank == right_rank: + mesh.comm.send(right_val, dest=0, tag=1) + + if mesh.comm.rank == 0: + left_val = mesh.comm.recv(source=left_rank, tag=0) + right_val = mesh.comm.recv(source=right_rank, tag=1) + assert np.isclose(J_sum, left_val + right_val) From b96f9acfa91f9336aea63605c5d4f40ee7f98d73 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= Date: Tue, 27 Aug 2024 12:46:00 +0000 Subject: [PATCH 2/4] Simplify MPI send recv --- python/test/unit/fem/test_assemble_submesh.py | 39 ++++++++++--------- 1 file changed, 20 insertions(+), 19 deletions(-) diff --git a/python/test/unit/fem/test_assemble_submesh.py b/python/test/unit/fem/test_assemble_submesh.py index 15e68bb042d..b16104b1920 100644 --- a/python/test/unit/fem/test_assemble_submesh.py +++ b/python/test/unit/fem/test_assemble_submesh.py @@ -523,31 +523,32 @@ def f_right(x): J_local = fem.assemble_scalar(J_compiled) J_sum = mesh.comm.allreduce(J_local, op=MPI.SUM) + vertex_map = mesh.topology.index_map(mesh.topology.dim-1) + num_vertices_local = vertex_map.size_local # Compute value of expression at left interface if len(facets := facet_tag.find(left_interface_tag)) > 0: + assert len(facets) == 1 left_vertex = entities_to_geometry(mesh, mesh.topology.dim - 1, facets) - left_coord = mesh.geometry.x[left_vertex].reshape(3, -1) - left_val = left_coord[0, 0] * f_left(left_coord)[0] + if left_vertex[0,0] < num_vertices_local: + left_coord = mesh.geometry.x[left_vertex].reshape(3, -1) + left_val = left_coord[0, 0] * f_left(left_coord)[0] + else: + left_val = 0.0 else: - left_val = 0 + left_val = 0.0 # Compute value of expression at right interface if len(facets := facet_tag.find(right_interface_tag)) > 0: + assert len(facets) == 1 right_vertex = entities_to_geometry(mesh, mesh.topology.dim - 1, facets) - right_coord = mesh.geometry.x[right_vertex].reshape(3, -1) - right_val = np.cos(right_coord[0, 0]) * f_right(right_coord)[0] + if right_vertex[0,0] < num_vertices_local: + right_coord = mesh.geometry.x[right_vertex].reshape(3, -1) + right_val = np.cos(right_coord[0, 0]) * f_right(right_coord)[0] + else: + right_val = 0.0 else: - right_val = 0 - - # As multiple process might have the vertex, only get value from one process - _, left_rank = mesh.comm.allreduce([abs(left_val), mesh.comm.rank], op=MPI.MAXLOC) - _, right_rank = mesh.comm.allreduce([abs(right_val), mesh.comm.rank], op=MPI.MAXLOC) - if mesh.comm.rank == left_rank: - mesh.comm.send(left_val, dest=0, tag=0) - if mesh.comm.rank == right_rank: - mesh.comm.send(right_val, dest=0, tag=1) - - if mesh.comm.rank == 0: - left_val = mesh.comm.recv(source=left_rank, tag=0) - right_val = mesh.comm.recv(source=right_rank, tag=1) - assert np.isclose(J_sum, left_val + right_val) + right_val = 0.0 + + glob_left_val = mesh.comm.allreduce(left_val,op=MPI.SUM) + glob_right_val = mesh.comm.allreduce(right_val,op=MPI.SUM) + assert np.isclose(J_sum, glob_left_val + glob_right_val) From f1fcc102f82dbed3ad2e56d9586b544953d37d24 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= Date: Tue, 27 Aug 2024 12:49:49 +0000 Subject: [PATCH 3/4] Ruff formatting --- python/test/unit/fem/test_assemble_submesh.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/python/test/unit/fem/test_assemble_submesh.py b/python/test/unit/fem/test_assemble_submesh.py index b16104b1920..ea768a2c8db 100644 --- a/python/test/unit/fem/test_assemble_submesh.py +++ b/python/test/unit/fem/test_assemble_submesh.py @@ -523,13 +523,13 @@ def f_right(x): J_local = fem.assemble_scalar(J_compiled) J_sum = mesh.comm.allreduce(J_local, op=MPI.SUM) - vertex_map = mesh.topology.index_map(mesh.topology.dim-1) - num_vertices_local = vertex_map.size_local + vertex_map = mesh.topology.index_map(mesh.topology.dim - 1) + num_vertices_local = vertex_map.size_local # Compute value of expression at left interface if len(facets := facet_tag.find(left_interface_tag)) > 0: assert len(facets) == 1 left_vertex = entities_to_geometry(mesh, mesh.topology.dim - 1, facets) - if left_vertex[0,0] < num_vertices_local: + if left_vertex[0, 0] < num_vertices_local: left_coord = mesh.geometry.x[left_vertex].reshape(3, -1) left_val = left_coord[0, 0] * f_left(left_coord)[0] else: @@ -541,7 +541,7 @@ def f_right(x): if len(facets := facet_tag.find(right_interface_tag)) > 0: assert len(facets) == 1 right_vertex = entities_to_geometry(mesh, mesh.topology.dim - 1, facets) - if right_vertex[0,0] < num_vertices_local: + if right_vertex[0, 0] < num_vertices_local: right_coord = mesh.geometry.x[right_vertex].reshape(3, -1) right_val = np.cos(right_coord[0, 0]) * f_right(right_coord)[0] else: @@ -549,6 +549,6 @@ def f_right(x): else: right_val = 0.0 - glob_left_val = mesh.comm.allreduce(left_val,op=MPI.SUM) - glob_right_val = mesh.comm.allreduce(right_val,op=MPI.SUM) + glob_left_val = mesh.comm.allreduce(left_val, op=MPI.SUM) + glob_right_val = mesh.comm.allreduce(right_val, op=MPI.SUM) assert np.isclose(J_sum, glob_left_val + glob_right_val) From 4834fb8dc349d3729fab0271ba3013aabe97a05e Mon Sep 17 00:00:00 2001 From: "Garth N. Wells" Date: Wed, 11 Sep 2024 10:49:26 +0100 Subject: [PATCH 4/4] Small logic simplification --- cpp/dolfinx/fem/utils.h | 46 ++++++++++--------- python/test/unit/fem/test_assemble_submesh.py | 13 +++--- 2 files changed, 32 insertions(+), 27 deletions(-) diff --git a/cpp/dolfinx/fem/utils.h b/cpp/dolfinx/fem/utils.h index c752e026fb4..cc1fda43e1e 100644 --- a/cpp/dolfinx/fem/utils.h +++ b/cpp/dolfinx/fem/utils.h @@ -843,7 +843,7 @@ get_cell_orientation_info(const Function& coefficient) } /// Pack a single coefficient for a single cell -template +template void pack(std::span coeffs, std::int32_t cell, int bs, std::span v, std::span cell_info, const DofMap& dofmap, auto transform) @@ -860,6 +860,7 @@ void pack(std::span coeffs, std::int32_t cell, int bs, std::span v, } else { + assert(_bs == bs); const int pos_c = _bs * i; const int pos_v = _bs * dofs[i]; for (int k = 0; k < _bs; ++k) @@ -917,44 +918,47 @@ void pack_coefficient_entity(std::span c, int cstride, for (std::size_t e = 0; e < entities.size(); e += estride) { auto entity = entities.subspan(e, estride); - std::int32_t cell = fetch_cells(entity); - if (cell < 0) - continue; - auto cell_coeff = c.subspan((e / estride) * cstride + offset, space_dim); - pack(cell_coeff, cell, bs, v, cell_info, dofmap, transformation); + if (std::int32_t cell = fetch_cells(entity); cell >= 0) + { + auto cell_coeff + = c.subspan((e / estride) * cstride + offset, space_dim); + pack<1>(cell_coeff, cell, bs, v, cell_info, dofmap, transformation); + } } break; case 2: for (std::size_t e = 0; e < entities.size(); e += estride) { auto entity = entities.subspan(e, estride); - std::int32_t cell = fetch_cells(entity); - if (cell < 0) - continue; - auto cell_coeff = c.subspan((e / estride) * cstride + offset, space_dim); - pack(cell_coeff, cell, bs, v, cell_info, dofmap, transformation); + if (std::int32_t cell = fetch_cells(entity); cell >= 0) + { + auto cell_coeff + = c.subspan((e / estride) * cstride + offset, space_dim); + pack<2>(cell_coeff, cell, bs, v, cell_info, dofmap, transformation); + } } break; case 3: for (std::size_t e = 0; e < entities.size(); e += estride) { auto entity = entities.subspan(e, estride); - std::int32_t cell = fetch_cells(entity); - if (cell < 0) - continue; - auto cell_coeff = c.subspan(e / estride * cstride + offset, space_dim); - pack(cell_coeff, cell, bs, v, cell_info, dofmap, transformation); + if (std::int32_t cell = fetch_cells(entity); cell >= 0) + { + auto cell_coeff = c.subspan(e / estride * cstride + offset, space_dim); + pack<3>(cell_coeff, cell, bs, v, cell_info, dofmap, transformation); + } } break; default: for (std::size_t e = 0; e < entities.size(); e += estride) { auto entity = entities.subspan(e, estride); - std::int32_t cell = fetch_cells(entity); - if (cell < 0) - continue; - auto cell_coeff = c.subspan((e / estride) * cstride + offset, space_dim); - pack(cell_coeff, cell, bs, v, cell_info, dofmap, transformation); + if (std::int32_t cell = fetch_cells(entity); cell >= 0) + { + auto cell_coeff + = c.subspan((e / estride) * cstride + offset, space_dim); + pack<-1>(cell_coeff, cell, bs, v, cell_info, dofmap, transformation); + } } break; } diff --git a/python/test/unit/fem/test_assemble_submesh.py b/python/test/unit/fem/test_assemble_submesh.py index ea768a2c8db..19056db2b30 100644 --- a/python/test/unit/fem/test_assemble_submesh.py +++ b/python/test/unit/fem/test_assemble_submesh.py @@ -442,17 +442,17 @@ def right(x): # One sided interface integral uses only "+" restriction. Sort integration entities such that # this is always satisfied def compute_mapped_interior_facet_data(mesh, facet_tag, value, parent_to_sub_map): - """ - Compute integration data for interior facet integrals, where the positive restriction is - always taken on the side that has a cell in the sub mesh + """Compute integration data for interior facet integrals, where the positive restriction is + always taken on the side that has a cell in the sub mesh. - Parameters + Args: mesh: Parent mesh facet_tag: Meshtags object for facets value: Value of the facets to extract parent_to_sub_map: Mapping from parent mesh to sub mesh - Returns - integration_data: Integration data for interior facets + + Returns: + Integration data for interior facets """ mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim) integration_data = compute_integration_domains( @@ -525,6 +525,7 @@ def f_right(x): vertex_map = mesh.topology.index_map(mesh.topology.dim - 1) num_vertices_local = vertex_map.size_local + # Compute value of expression at left interface if len(facets := facet_tag.find(left_interface_tag)) > 0: assert len(facets) == 1