From 78dc73df9dd954b680d6a1d4ab6f6519049194e5 Mon Sep 17 00:00:00 2001 From: Reuben Nixon-Hill Date: Fri, 20 Jan 2023 13:36:52 +0000 Subject: [PATCH 01/14] Correctly name newton step function to_reference_coordinates is an inaccurate name for the function, it actually only does a single newton step! --- firedrake/mg/kernels.py | 6 +++--- firedrake/pointquery_utils.py | 6 +++--- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/firedrake/mg/kernels.py b/firedrake/mg/kernels.py index 3d2e2d7da0..6984f26d50 100644 --- a/firedrake/mg/kernels.py +++ b/firedrake/mg/kernels.py @@ -29,7 +29,7 @@ from tsfc.finatinterface import create_element from finat.quadrature import make_quadrature from firedrake.pointquery_utils import dX_norm_square, X_isub_dX, init_X, inside_check, is_affine, compute_celldist -from firedrake.pointquery_utils import to_reference_coordinates as to_reference_coordinates_body +from firedrake.pointquery_utils import to_reference_coords_newton_step as to_reference_coords_newton_step_body def to_reference_coordinates(ufl_coordinate_element, parameters=None): @@ -48,7 +48,7 @@ def to_reference_coordinates(ufl_coordinate_element, parameters=None): code = { "geometric_dimension": cell.geometric_dimension(), "topological_dimension": cell.topological_dimension(), - "to_reference_coords": to_reference_coordinates_body(ufl_coordinate_element, parameters), + "to_reference_coords_newton_step": to_reference_coords_newton_step_body(ufl_coordinate_element, parameters), "init_X": init_X(element.cell, parameters), "max_iteration_count": 1 if is_affine(ufl_coordinate_element) else 16, "convergence_epsilon": 1e-12, @@ -75,7 +75,7 @@ def to_reference_coordinates(ufl_coordinate_element, parameters=None): int converged = 0; for (int it = 0; !converged && it < %(max_iteration_count)d; it++) { double dX[%(topological_dimension)d] = { 0.0 }; -%(to_reference_coords)s +%(to_reference_coords_newton_step)s if (%(dX_norm_square)s < %(convergence_epsilon)g * %(convergence_epsilon)g) { converged = 1; diff --git a/firedrake/pointquery_utils.py b/firedrake/pointquery_utils.py index 73bae09b52..11f26c4bd6 100644 --- a/firedrake/pointquery_utils.py +++ b/firedrake/pointquery_utils.py @@ -95,7 +95,7 @@ def init_X(fiat_cell, parameters): @PETSc.Log.EventDecorator() -def to_reference_coordinates(ufl_coordinate_element, parameters): +def to_reference_coords_newton_step(ufl_coordinate_element, parameters): # Set up UFL form cell = ufl_coordinate_element.cell() domain = ufl.Mesh(ufl_coordinate_element) @@ -172,7 +172,7 @@ def compile_coordinate_element(ufl_coordinate_element, contains_eps, parameters= "geometric_dimension": cell.geometric_dimension(), "topological_dimension": cell.topological_dimension(), "inside_predicate": inside_check(element.cell, eps=contains_eps), - "to_reference_coords": to_reference_coordinates(ufl_coordinate_element, parameters), + "to_reference_coords_newton_step": to_reference_coords_newton_step(ufl_coordinate_element, parameters), "init_X": init_X(element.cell, parameters), "max_iteration_count": 1 if is_affine(ufl_coordinate_element) else 16, "convergence_epsilon": 1e-12, @@ -204,7 +204,7 @@ def compile_coordinate_element(ufl_coordinate_element, contains_eps, parameters= int converged = 0; for (int it = 0; !converged && it < %(max_iteration_count)d; it++) { %(ScalarType)s dX[%(topological_dimension)d] = { 0.0 }; -%(to_reference_coords)s +%(to_reference_coords_newton_step)s if (%(dX_norm_square)s < %(convergence_epsilon)g * %(convergence_epsilon)g) { converged = 1; From 588e35cd8d264ab030e4bc8967e792c44b6b39c1 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Mon, 23 Jan 2023 14:27:00 +0000 Subject: [PATCH 02/14] Explicitly check whether parmetis is installed as a PETSc package --- firedrake_configuration/__init__.py | 15 ++++++++++++++- 1 file changed, 14 insertions(+), 1 deletion(-) diff --git a/firedrake_configuration/__init__.py b/firedrake_configuration/__init__.py index f6e49857d4..b7970ab0b7 100644 --- a/firedrake_configuration/__init__.py +++ b/firedrake_configuration/__init__.py @@ -7,6 +7,7 @@ import json import os import sys +import petsc4py # Attempt to read configuration from file. try: @@ -22,7 +23,19 @@ _config = json.load(f) except IOError: - _config = None + _config = {} + + +def petsc_packages(): + conf = petsc4py.get_config() + with open(os.path.join(conf["PETSC_DIR"], conf["PETSC_ARCH"], "include", "petscconf.h"), "r") as f: + *_, packages = next(line for line in f if line.startswith("#define PETSC_HAVE_PACKAGES")).split() + return set(packages[2:-2].split(":")) + + +options = _config.get("options", {}) +options["with_parmetis"] = "parmetis" in petsc_packages() +_config["options"] = options def get_config(): From 6115e2f2c4d1f8654a0bd46257af4da5a127b16f Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Mon, 23 Jan 2023 15:24:06 +0000 Subject: [PATCH 03/14] Fix DirichletBCs on diagonal of mixed bilinear --- firedrake/bcs.py | 5 ++++- tests/regression/test_bcs.py | 15 +++++++++------ 2 files changed, 13 insertions(+), 7 deletions(-) diff --git a/firedrake/bcs.py b/firedrake/bcs.py index a7c01295b0..c660585c5a 100644 --- a/firedrake/bcs.py +++ b/firedrake/bcs.py @@ -193,9 +193,12 @@ def set(self, r, val): :arg r: the :class:`Function` to which the value should be applied. :arg val: the prescribed value. """ + for idx in self._indices: r = r.sub(idx) - val = val.sub(idx) + if not np.isscalar(val): + for idx in self._indices: + val = val.sub(idx) r.assign(val, subset=self.node_set) def integrals(self): diff --git a/tests/regression/test_bcs.py b/tests/regression/test_bcs.py index deef1e5015..d55404a4a8 100644 --- a/tests/regression/test_bcs.py +++ b/tests/regression/test_bcs.py @@ -276,7 +276,9 @@ def test_overlapping_bc_nodes(quad): assert np.allclose(A, np.identity(V.dof_dset.size)) -def test_mixed_bcs(): +@pytest.mark.parametrize("diagonal", + [False, True]) +def test_mixed_bcs(diagonal): m = UnitSquareMesh(2, 2) V = FunctionSpace(m, 'CG', 1) W = V*V @@ -285,11 +287,12 @@ def test_mixed_bcs(): bc = DirichletBC(W.sub(1), 0.0, "on_boundary") - A = assemble(inner(u, v)*dx, bcs=bc) - - A11 = A.M[1, 1].values - - assert np.allclose(A11.diagonal()[bc.nodes], 1.0) + A = assemble(inner(u, v)*dx, bcs=bc, diagonal=diagonal) + if diagonal: + data = A.dat[1].data + else: + data = A.M[1, 1].values.diagonal() + assert np.allclose(data[bc.nodes], 1.0) def test_bcs_rhs_assemble(a, V): From a149acb362744f52745ae0a2fe84f432c4b3893b Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Tue, 24 Jan 2023 09:28:38 +0000 Subject: [PATCH 04/14] Set ids on test --- tests/regression/test_bcs.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/tests/regression/test_bcs.py b/tests/regression/test_bcs.py index d55404a4a8..24c9d56b46 100644 --- a/tests/regression/test_bcs.py +++ b/tests/regression/test_bcs.py @@ -277,7 +277,8 @@ def test_overlapping_bc_nodes(quad): @pytest.mark.parametrize("diagonal", - [False, True]) + [False, True], + ids=["matrix", "diagonal"]) def test_mixed_bcs(diagonal): m = UnitSquareMesh(2, 2) V = FunctionSpace(m, 'CG', 1) From ad4d08025c43f8b160e2a771f6a194181d696a4a Mon Sep 17 00:00:00 2001 From: Jack Betteridge <43041811+JDBetteridge@users.noreply.github.com> Date: Wed, 25 Jan 2023 11:27:38 +0000 Subject: [PATCH 05/14] Cancel build on new commit (#2723) --- .github/workflows/build.yml | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml index 6287cb38d6..cacb353799 100644 --- a/.github/workflows/build.yml +++ b/.github/workflows/build.yml @@ -9,6 +9,11 @@ on: schedule: - cron: '0 0 * * 0' +concurrency: + # Cancels jobs running if new commits are pushed + group: ${{ github.workflow }}-${{ github.event.pull_request.number || github.ref }} + cancel-in-progress: true + jobs: build: name: "Build Firedrake" From 5fc31e85737a13bdae5cc16151be5aefe3da0c2c Mon Sep 17 00:00:00 2001 From: Jack Betteridge <43041811+JDBetteridge@users.noreply.github.com> Date: Wed, 25 Jan 2023 11:28:28 +0000 Subject: [PATCH 06/14] Fix tigerlake gcc 9 detection (#2726) * Fix tigerlake gcc 9 detection * Skip tigerlake gcc 9 detection on update * Tidy code --- scripts/firedrake-install | 28 +++++++++++++++------------- 1 file changed, 15 insertions(+), 13 deletions(-) diff --git a/scripts/firedrake-install b/scripts/firedrake-install index a9f78d49e4..afb39f5baf 100755 --- a/scripts/firedrake-install +++ b/scripts/firedrake-install @@ -518,23 +518,25 @@ def sniff_compiler(exe): else: name = "unknown" # VERSION - try: - output = subprocess.run( - [exe, "-dumpversion"], - stdout=subprocess.PIPE, - stderr=subprocess.PIPE, - check=True, - encoding="utf-8" - ).stdout - version = Version(output) - except (subprocess.CalledProcessError, UnicodeDecodeError, InvalidVersion): - version = None + version = None + for dumpstring in ["-dumpfullversion", "-dumpversion"]: + try: + output = subprocess.run( + [exe, dumpstring], + stdout=subprocess.PIPE, + stderr=subprocess.PIPE, + check=True, + encoding="utf-8" + ).stdout + version = Version(output) + break + except (subprocess.CalledProcessError, UnicodeDecodeError, InvalidVersion): + continue return name, version # Temporary catch for tigerlake processors with gcc/g++/gfortran 9.3 -# But only for Linux -if platform.system() == "Linux": +if platform.system() == "Linux" and mode == "install": # Get the cpu information try: cpuinfostring = subprocess.run( From 543e9d30f2f9b215857e6a2c9dc16d878fd46efe Mon Sep 17 00:00:00 2001 From: Nacime Bouziani <48448063+nbouziani@users.noreply.github.com> Date: Wed, 25 Jan 2023 11:37:55 +0000 Subject: [PATCH 07/14] Check UFL FEniCS merge (#2715) * Update build.yml * Remove deprecated ufl elements * Update tests * Fix multiple subdomain data consequences * Update few files using deprecated features * Restore and update tests for facet and interior elements * Clean notebook diff * Update build.yml --- .../higher_order_mass_lumping.py.rst | 2 +- docs/notebooks/10-sum-factorisation.ipynb | 6 ++-- firedrake/assemble.py | 28 +++++++++++-------- firedrake/functionspace.py | 3 +- firedrake/slate/slate.py | 11 ++++---- firedrake/utility_meshes.py | 2 +- tests/regression/test_facet_elements.py | 2 +- tests/regression/test_function_spaces.py | 3 +- tests/regression/test_hyperelasticity.py | 2 +- tests/regression/test_interior_elements.py | 2 +- tests/regression/test_project_interp_KMV.py | 2 +- tests/slate/test_cg_poisson.py | 2 +- tests/slate/test_slate_infrastructure.py | 27 ++++++++++++------ 13 files changed, 54 insertions(+), 38 deletions(-) diff --git a/demos/higher_order_mass_lumping/higher_order_mass_lumping.py.rst b/demos/higher_order_mass_lumping/higher_order_mass_lumping.py.rst index c22cf8e2dd..aca4f11032 100644 --- a/demos/higher_order_mass_lumping/higher_order_mass_lumping.py.rst +++ b/demos/higher_order_mass_lumping/higher_order_mass_lumping.py.rst @@ -114,7 +114,7 @@ space, along with the degree of the element and construct the quadrature rule:: Then we make a new Measure object that uses this rule:: - dxlump=dx(rule=quad_rule) + dxlump=dx(scheme=quad_rule) To discretize :math:`\partial_{t}^2 u` we use a central scheme diff --git a/docs/notebooks/10-sum-factorisation.ipynb b/docs/notebooks/10-sum-factorisation.ipynb index 63d07ec78f..4f942ebff6 100644 --- a/docs/notebooks/10-sum-factorisation.ipynb +++ b/docs/notebooks/10-sum-factorisation.ipynb @@ -265,7 +265,7 @@ "outputs": [], "source": [ "gll_quadrature_rule = gauss_lobatto_legendre_cube_rule(dimension=3, degree=p)\n", - "a_gll = dot(grad(u), grad(v)) *dx(rule=gll_quadrature_rule)" + "a_gll = dot(grad(u), grad(v)) *dx(scheme=gll_quadrature_rule)" ] }, { @@ -326,7 +326,7 @@ " V = FunctionSpace(mesh, element)\n", " u = TrialFunction(V)\n", " v = TestFunction(V)\n", - " a = dot(grad(u), grad(v))*dx(rule=modes[mode]['rule'](3, p))\n", + " a = dot(grad(u), grad(v))*dx(scheme=modes[mode]['rule'](3, p))\n", " kernel, = compile_form(a, parameters={\"mode\": modes[mode]['mode']})\n", " flops[mode].append(EstimateFlops().visit(kernel.ast))" ] @@ -1163,7 +1163,7 @@ " V = FunctionSpace(mesh, element)\n", " u = TrialFunction(V)\n", " v = TestFunction(V)\n", - " a = dot(curl(u), curl(v))*dx(rule=modes[mode]['rule'](3, p))\n", + " a = dot(curl(u), curl(v))*dx(scheme=modes[mode]['rule'](3, p))\n", " kernel, = compile_form(a, parameters={\"mode\": modes[mode]['mode']})\n", " flops_curl[mode].append(EstimateFlops().visit(kernel.ast))" ] diff --git a/firedrake/assemble.py b/firedrake/assemble.py index 9c9ab6dc24..de1aec7792 100644 --- a/firedrake/assemble.py +++ b/firedrake/assemble.py @@ -730,13 +730,14 @@ def _global_kernel_cache_key(form, local_knl, all_integer_subdomain_ids, **kwarg subdomain_key = [] for val in form.subdomain_data().values(): for k, v in val.items(): - if v is not None: - extruded = v._extruded - constant_layers = extruded and v.constant_layers - subset = isinstance(v, op2.Subset) - subdomain_key.append((k, extruded, constant_layers, subset)) - else: - subdomain_key.append((k,)) + for i, vi in enumerate(v): + if vi is not None: + extruded = vi._extruded + constant_layers = extruded and vi.constant_layers + subset = isinstance(vi, op2.Subset) + subdomain_key.append((k, i, extruded, constant_layers, subset)) + else: + subdomain_key.append((k, i)) return ((sig,) + tuple(subdomain_key) @@ -812,7 +813,7 @@ def _mesh(self): @cached_property def _needs_subset(self): subdomain_data = self._form.subdomain_data()[self._mesh] - if subdomain_data.get(self._integral_type, None) is not None: + if not all(sd is None for sd in subdomain_data.get(self._integral_type, [None])): return True if self._kinfo.subdomain_id == "everywhere": @@ -1072,9 +1073,14 @@ def _iterset(self): try: subdomain_data = self._form.subdomain_data()[self._mesh][self._integral_type] except KeyError: - subdomain_data = None - - if subdomain_data is not None: + subdomain_data = [None] + + subdomain_data = [sd for sd in subdomain_data if sd is not None] + if subdomain_data: + try: + subdomain_data, = subdomain_data + except ValueError: + raise NotImplementedError("Assembly with multiple subdomain data values is not supported") if self._integral_type != "cell": raise NotImplementedError("subdomain_data only supported with cell integrals") if self._kinfo.subdomain_id not in ["everywhere", "otherwise"]: diff --git a/firedrake/functionspace.py b/firedrake/functionspace.py index fc7592acca..01f154bcb1 100644 --- a/firedrake/functionspace.py +++ b/firedrake/functionspace.py @@ -78,8 +78,7 @@ def check_element(element, top=True): :raises ValueError: if the element is illegal. """ - if type(element) in (ufl.BrokenElement, ufl.FacetElement, - ufl.InteriorElement, ufl.RestrictedElement, + if type(element) in (ufl.BrokenElement, ufl.RestrictedElement, ufl.HDivElement, ufl.HCurlElement): inner = (element._element, ) elif type(element) is ufl.EnrichedElement: diff --git a/firedrake/slate/slate.py b/firedrake/slate/slate.py index 169ec0c07d..c70bc35870 100644 --- a/firedrake/slate/slate.py +++ b/firedrake/slate/slate.py @@ -478,7 +478,7 @@ def subdomain_data(self): """Returns a mapping on the tensor: ``{domain:{integral_type: subdomain_data}}``. """ - return {self.ufl_domain(): {"cell": None}} + return {self.ufl_domain(): {"cell": [None]}} def _output_string(self, prec=None): """Creates a string representation of the tensor.""" @@ -553,7 +553,7 @@ def subdomain_data(self): """Returns mappings on the tensor: ``{domain:{integral_type: subdomain_data}}``. """ - return tuple({domain: {"cell": None}} for domain in self.ufl_domain()) + return tuple({domain: {"cell": [None]}} for domain in self.ufl_domain()) def _output_string(self, prec=None): """Creates a string representation of the tensor.""" @@ -970,9 +970,10 @@ def subdomain_data(self): sd[it_type] = domain else: - assert sd[it_type] == domain, ( - "Domains must agree!" - ) + if not all(d is None for d in sd[it_type]) or not all(d is None for d in domain): + assert sd[it_type] == domain, ( + "Domains must agree!" + ) return {self.ufl_domain(): sd} diff --git a/firedrake/utility_meshes.py b/firedrake/utility_meshes.py index 9e39e1d43a..1626c036e1 100644 --- a/firedrake/utility_meshes.py +++ b/firedrake/utility_meshes.py @@ -2112,7 +2112,7 @@ def OctahedralSphereMesh( # This will DTWT on meshes with more than 26 refinement levels. # (log_2 1e8 ~= 26.5) tol = ufl.real(Constant(1.0e-8)) - rnew = ufl.Max(1 - abs(z), 0) + rnew = ufl.max_value(1 - abs(z), 0) # Avoid division by zero (when rnew is zero, x & y are also zero) x0 = ufl.conditional(ufl.lt(ufl.real(rnew), tol), 0, x / rnew) y0 = ufl.conditional(ufl.lt(rnew, tol), 0, y / rnew) diff --git a/tests/regression/test_facet_elements.py b/tests/regression/test_facet_elements.py index 46a3e0beeb..6e368dd68c 100644 --- a/tests/regression/test_facet_elements.py +++ b/tests/regression/test_facet_elements.py @@ -5,7 +5,7 @@ def test_all_dofs_on_facets(): mesh = UnitSquareMesh(5, 5) V_elt = FiniteElement("BDM", triangle, 1) # BDM has all dofs on facets, so these should be the same - W1_elt = FacetElement(V_elt) + W1_elt = RestrictedElement(V_elt, "facet") V = FunctionSpace(mesh, V_elt) W1 = FunctionSpace(mesh, W1_elt) f = Function(V).assign(1) diff --git a/tests/regression/test_function_spaces.py b/tests/regression/test_function_spaces.py index 73ddfa2600..1d33c0f0d6 100644 --- a/tests/regression/test_function_spaces.py +++ b/tests/regression/test_function_spaces.py @@ -95,8 +95,7 @@ def test_function_space_collapse(cg1): @pytest.mark.parametrize("modifier", - [BrokenElement, FacetElement, - InteriorElement, HDivElement, + [BrokenElement, HDivElement, HCurlElement]) @pytest.mark.parametrize("element", [FiniteElement("CG", triangle, 1), diff --git a/tests/regression/test_hyperelasticity.py b/tests/regression/test_hyperelasticity.py index eac7c0e0a2..a2b7125057 100644 --- a/tests/regression/test_hyperelasticity.py +++ b/tests/regression/test_hyperelasticity.py @@ -27,7 +27,7 @@ def test_hyperelastic_convergence(): T = Constant((0, -0.5)) B = Constant((0, -0.25)) - d = u.geometric_dimension() + d = mesh.geometric_dimension() I = Identity(d) F = I + grad(u) # Deformation gradient C = F.T*F # Right Cauchy-Green tensor diff --git a/tests/regression/test_interior_elements.py b/tests/regression/test_interior_elements.py index 885a7848b9..faf880cbae 100644 --- a/tests/regression/test_interior_elements.py +++ b/tests/regression/test_interior_elements.py @@ -4,7 +4,7 @@ def test_vanish_on_bdy(): mesh = UnitSquareMesh(5, 5) V_elt = FiniteElement("RT", triangle, 2) - W2_elt = InteriorElement(V_elt) + W2_elt = RestrictedElement(V_elt, "interior") W2 = FunctionSpace(mesh, W2_elt) g = Function(W2).assign(1) n = FacetNormal(mesh) diff --git a/tests/regression/test_project_interp_KMV.py b/tests/regression/test_project_interp_KMV.py index 636b6fefe1..f3b0ae600d 100644 --- a/tests/regression/test_project_interp_KMV.py +++ b/tests/regression/test_project_interp_KMV.py @@ -60,7 +60,7 @@ def run_projection(mesh, expr, p): r = Function(V) f = interpolate(expr(*SpatialCoordinate(mesh)), V) solve( - inner(u, v) * dx(rule=qr) == inner(f, v) * dx(rule=qr), + inner(u, v) * dx(scheme=qr) == inner(f, v) * dx(scheme=qr), r, solver_parameters={ "ksp_type": "preonly", diff --git a/tests/slate/test_cg_poisson.py b/tests/slate/test_cg_poisson.py index f04d545ef1..4d9665914a 100644 --- a/tests/slate/test_cg_poisson.py +++ b/tests/slate/test_cg_poisson.py @@ -25,7 +25,7 @@ def run_CG_problem(r, degree, quads=False): # Set up function spaces e = FiniteElement("Lagrange", cell=mesh.ufl_cell(), degree=degree) - Z = FunctionSpace(mesh, MixedElement(InteriorElement(e), FacetElement(e))) + Z = FunctionSpace(mesh, MixedElement(RestrictedElement(e, "interior"), RestrictedElement(e, "facet"))) z = Function(Z) u = sum(split(z)) diff --git a/tests/slate/test_slate_infrastructure.py b/tests/slate/test_slate_infrastructure.py index db3aa33ab5..156bb3caff 100644 --- a/tests/slate/test_slate_infrastructure.py +++ b/tests/slate/test_slate_infrastructure.py @@ -195,14 +195,14 @@ def test_integral_information(mass, stiffness, load, boundary_load, zero_rank_te assert (F + G).ufl_domain() == (F.form + G.form).ufl_domain() assert (M + N).ufl_domain() == (M.form + N.form).ufl_domain() - assert S.subdomain_data() == S.form.subdomain_data() - assert N.subdomain_data() == N.form.subdomain_data() - assert M.subdomain_data() == M.form.subdomain_data() - assert F.subdomain_data() == F.form.subdomain_data() - assert N.inv.subdomain_data() == N.form.subdomain_data() - assert (-M).subdomain_data() == M.form.subdomain_data() - assert (M + N).T.subdomain_data() == (M.form + N.form).subdomain_data() - assert (F + G).subdomain_data() == (F.form + G.form).subdomain_data() + assert _is_equal_subdomain_data(S, S.form) + assert _is_equal_subdomain_data(N, N.form) + assert _is_equal_subdomain_data(M, M.form) + assert _is_equal_subdomain_data(F, F.form) + assert _is_equal_subdomain_data(N.inv, N.form) + assert _is_equal_subdomain_data(-M, M.form) + assert _is_equal_subdomain_data((M + N).T, (M.form + N.form)) + assert _is_equal_subdomain_data((F + G), (F.form + G.form)) def test_equality_relations(function_space): @@ -387,3 +387,14 @@ def test_illegal_compile(): form = v * dx with pytest.raises(ValueError): compile_slate(form) + + +def _is_equal_subdomain_data(a, b): + """Compare subdomain data of a and b.""" + sd_a = {domain: {integral_type: [v for v in val if v is not None]} + for domain, data in a.subdomain_data().items() + for integral_type, val in data.items()} + sd_b = {domain: {integral_type: [v for v in val if v is not None]} + for domain, data in b.subdomain_data().items() + for integral_type, val in data.items()} + return sd_a == sd_b From bc5fa32b059cd640789cef4c97b430e4ea473824 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Wed, 25 Jan 2023 16:11:32 +0000 Subject: [PATCH 08/14] Projector: matfree solver params (#2728) * Projection: set Jacobi as default solver * test mixed projector with aij, nest, and matfree --- firedrake/projection.py | 12 +++++++-- tests/regression/test_projection.py | 39 +++++++++++++++++++++++++++-- 2 files changed, 47 insertions(+), 4 deletions(-) diff --git a/firedrake/projection.py b/firedrake/projection.py index 05c697400c..475646a310 100644 --- a/firedrake/projection.py +++ b/firedrake/projection.py @@ -99,8 +99,16 @@ def __init__(self, source, target, bcs=None, solver_parameters=None, solver_parameters = solver_parameters.copy() solver_parameters.setdefault("ksp_type", "cg") solver_parameters.setdefault("ksp_rtol", 1e-8) - solver_parameters.setdefault("pc_type", "bjacobi") - solver_parameters.setdefault("sub_pc_type", "icc") + mat_type = solver_parameters.get("mat_type", firedrake.parameters["default_matrix_type"]) + if mat_type == "nest": + solver_parameters.setdefault("pc_type", "fieldsplit") + solver_parameters.setdefault("fieldsplit_pc_type", "bjacobi") + solver_parameters.setdefault("fieldsplit_sub_pc_type", "icc") + elif mat_type == "matfree": + solver_parameters.setdefault("pc_type", "jacobi") + else: + solver_parameters.setdefault("pc_type", "bjacobi") + solver_parameters.setdefault("sub_pc_type", "icc") self.source = source self.target = target self.solver_parameters = solver_parameters diff --git a/tests/regression/test_projection.py b/tests/regression/test_projection.py index b8d28b1e98..8efe9f453a 100644 --- a/tests/regression/test_projection.py +++ b/tests/regression/test_projection.py @@ -193,7 +193,8 @@ def test_repeatable(): assert (fd == ud).all() -def test_projector(): +@pytest.mark.parametrize('mat_type', ['aij', 'matfree']) +def test_projector(mat_type): m = UnitSquareMesh(2, 2) Vc = FunctionSpace(m, "CG", 2) xs = SpatialCoordinate(m) @@ -204,7 +205,7 @@ def test_projector(): Vd = FunctionSpace(m, "DG", 1) vo = Function(Vd) - P = Projector(v, vo) + P = Projector(v, vo, solver_parameters={"mat_type": mat_type}) P.project() mass2 = assemble(vo*dx) @@ -218,6 +219,40 @@ def test_projector(): assert np.abs(mass1-mass2) < 1.0e-10 +@pytest.mark.parametrize('mat_type', ['aij', 'nest', 'matfree']) +def test_mixed_projector(mat_type): + m = UnitSquareMesh(2, 2) + Vc1 = FunctionSpace(m, "CG", 1) + Vc2 = FunctionSpace(m, "CG", 2) + Vc = Vc1 * Vc2 + xs = SpatialCoordinate(m) + + v = Function(Vc) + v0, v1 = v.split() # sub_functions + v0.interpolate(xs[0]*xs[1] + cos(xs[0]+xs[1])) + v1.interpolate(xs[0]*xs[1] + sin(xs[0]+xs[1])) + mass1 = assemble(sum(split(v))*dx) + + Vd1 = FunctionSpace(m, "DG", 1) + Vd2 = FunctionSpace(m, "DG", 2) + Vd = Vd1 * Vd2 + vo = Function(Vd) + + P = Projector(v, vo, solver_parameters={"mat_type": mat_type}) + P.project() + + mass2 = assemble(sum(split(vo))*dx) + assert np.abs(mass1-mass2) < 1.0e-10 + + v0.interpolate(xs[1] + exp(xs[0]+xs[1])) + v1.interpolate(xs[0] + exp(xs[0]+xs[1])) + mass1 = assemble(sum(split(v))*dx) + + P.project() + mass2 = assemble(sum(split(vo))*dx) + assert np.abs(mass1-mass2) < 1.0e-10 + + def test_trivial_projector(): m = UnitSquareMesh(2, 2) Vc = FunctionSpace(m, "CG", 2) From 01d111ab3818904df318b4441fbaa0cd9dd80e90 Mon Sep 17 00:00:00 2001 From: Reuben Nixon-Hill Date: Thu, 26 Jan 2023 11:15:17 +0000 Subject: [PATCH 09/14] use sympys in built ccode generation --- firedrake/pointquery_utils.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/firedrake/pointquery_utils.py b/firedrake/pointquery_utils.py index 11f26c4bd6..da9f0a4fdc 100644 --- a/firedrake/pointquery_utils.py +++ b/firedrake/pointquery_utils.py @@ -1,6 +1,7 @@ from os import path import numpy import sympy +from sympy.printing.c import ccode from pyop2 import op2 from pyop2.parloop import generate_single_cell_wrapper @@ -65,8 +66,7 @@ def is_affine(ufl_element): def inside_check(fiat_cell, eps, X="X"): dim = fiat_cell.get_spatial_dimension() point = tuple(sympy.Symbol("PetscRealPart(%s[%d])" % (X, i)) for i in range(dim)) - - return " && ".join("(%s)" % arg for arg in fiat_cell.contains_point(point, epsilon=eps).args) + return ccode(fiat_cell.contains_point(point, epsilon=eps)) def compute_celldist(fiat_cell, X="X", celldist="celldist"): From e63871794dae2516bff6820d5aeaf544e863d70f Mon Sep 17 00:00:00 2001 From: Reuben Nixon-Hill Date: Mon, 30 Jan 2023 11:24:53 +0000 Subject: [PATCH 10/14] add missing tolerance documentation --- firedrake/mesh.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/firedrake/mesh.py b/firedrake/mesh.py index 9b71321fbc..7a388e7d99 100644 --- a/firedrake/mesh.py +++ b/firedrake/mesh.py @@ -2435,6 +2435,10 @@ def _pic_swarm_in_mesh(parent_mesh, coords, fields=None, tolerance=None): RealType)]``. All fields must have the same number of points. For more information see `the DMSWARM API reference _. + :kwarg tolerance: The tolerance used by locate_cell when deciding which + cell each DMSwarm point is found. The default is None which can cause + problems if the DMSwarm points are at the boundary of two cells or just + outside the mesh. :return: the immersed DMSwarm .. note:: From f9acbc048e91bf193db7ab2ef6b8edf459ac2cb1 Mon Sep 17 00:00:00 2001 From: Reuben Nixon-Hill Date: Fri, 27 Jan 2023 09:57:25 +0000 Subject: [PATCH 11/14] have a default tolerance to match docstring --- firedrake/mesh.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/firedrake/mesh.py b/firedrake/mesh.py index 7a388e7d99..d39a1ed79d 100644 --- a/firedrake/mesh.py +++ b/firedrake/mesh.py @@ -2250,7 +2250,7 @@ def ExtrudedMesh(mesh, layers, layer_height=None, extrusion_type='uniform', kern @PETSc.Log.EventDecorator() def VertexOnlyMesh(mesh, vertexcoords, missing_points_behaviour=None, - tolerance=None): + tolerance=1.0e-14): """ Create a vertex only mesh, immersed in a given mesh, with vertices defined by a list of coordinates. From b79c2d0d33e53b7e7f4a28eb05762c6c14018874 Mon Sep 17 00:00:00 2001 From: Reuben Nixon-Hill Date: Fri, 27 Jan 2023 17:25:41 +0000 Subject: [PATCH 12/14] Add default redundant option to VertexOnlyMesh This causes the mesh to be built using the points supplied on rank 0 only and is switched on by default. Various tests are also present. --- firedrake/mesh.py | 32 +++++++++------ tests/vertexonly/test_swarm.py | 39 +++++++++++++++---- .../test_vertex_only_mesh_generation.py | 31 +++++++++++++-- 3 files changed, 80 insertions(+), 22 deletions(-) diff --git a/firedrake/mesh.py b/firedrake/mesh.py index d39a1ed79d..4a8334865e 100644 --- a/firedrake/mesh.py +++ b/firedrake/mesh.py @@ -2250,7 +2250,7 @@ def ExtrudedMesh(mesh, layers, layer_height=None, extrusion_type='uniform', kern @PETSc.Log.EventDecorator() def VertexOnlyMesh(mesh, vertexcoords, missing_points_behaviour=None, - tolerance=1.0e-14): + tolerance=1.0e-14, redundant=True): """ Create a vertex only mesh, immersed in a given mesh, with vertices defined by a list of coordinates. @@ -2269,6 +2269,11 @@ def VertexOnlyMesh(mesh, vertexcoords, missing_points_behaviour=None, in the cell. Increase the default (1.0e-14) somewhat if vertices interior to the domain are being lost in the :class:`VertexOnlyMesh` construction process. + :kwarg redundant: If True, the mesh will be built using just the vertices + which are specified on rank 0. If False, the mesh will be built using + the vertices specified by each rank. Care must be taken when using + ``redundant = False``: see the note below for more information. + .. note:: @@ -2285,17 +2290,17 @@ def VertexOnlyMesh(mesh, vertexcoords, missing_points_behaviour=None, the VertexOnlyMesh to return the wrong values. .. note:: - When running in parallel, ``vertexcoords`` are strictly confined - to the local ``mesh`` cells of that rank. This means that if rank - A has ``vertexcoords`` {X} that are not found in the mesh cells - owned by rank A but are found in the mesh cells owned by rank B, - **and rank B has not been supplied with those** ``vertexcoords``, - then the ``vertexcoords`` {X} will be lost. + When running in parallel with ``redundant = False``, ``vertexcoords`` + are strictly confined to the local ``mesh`` cells of that rank. This + means that if rank A has ``vertexcoords`` {X} that are not found in the + mesh cells owned by rank A but are found in the mesh cells owned by + ank B, **and rank B has not been supplied with those** + ``vertexcoords``, then the ``vertexcoords`` {X} will be lost. This can be avoided by either - #. making sure that all ranks are supplied with the same - ``vertexcoords`` or by + #. making sure that rank 0 has all vertex coordinates and setting + ``redundant = True`` or by #. ensuring that ``vertexcoords`` are already found in cells owned by the ``mesh`` partition of the given rank. @@ -2339,7 +2344,7 @@ def VertexOnlyMesh(mesh, vertexcoords, missing_points_behaviour=None, if pdim != gdim: raise ValueError(f"Mesh geometric dimension {gdim} must match point list dimension {pdim}") - swarm = _pic_swarm_in_mesh(mesh, vertexcoords, tolerance=tolerance) + swarm = _pic_swarm_in_mesh(mesh, vertexcoords, tolerance=tolerance, redundant=redundant) if missing_points_behaviour: @@ -2413,7 +2418,7 @@ def compare_arrays(x, y, datatype): return vmesh -def _pic_swarm_in_mesh(parent_mesh, coords, fields=None, tolerance=None): +def _pic_swarm_in_mesh(parent_mesh, coords, fields=None, tolerance=None, redundant=True): """Create a Particle In Cell (PIC) DMSwarm immersed in a Mesh This should only by used for meshes with straight edges. If not, the @@ -2439,6 +2444,8 @@ def _pic_swarm_in_mesh(parent_mesh, coords, fields=None, tolerance=None): cell each DMSwarm point is found. The default is None which can cause problems if the DMSwarm points are at the boundary of two cells or just outside the mesh. + :kwarg redundant: If True, the DMSwarm will be created using only the + points specified on MPI rank 0. The default is True. :return: the immersed DMSwarm .. note:: @@ -2473,6 +2480,9 @@ def _pic_swarm_in_mesh(parent_mesh, coords, fields=None, tolerance=None): # Check coords coords = np.asarray(coords, dtype=RealType) + if redundant: + coords = parent_mesh._comm.bcast(coords, root=0) + plex = parent_mesh.topology.topology_dm tdim = parent_mesh.topological_dimension() gdim = parent_mesh.geometric_dimension() diff --git a/tests/vertexonly/test_swarm.py b/tests/vertexonly/test_swarm.py index 152eef3b55..2f48875b80 100644 --- a/tests/vertexonly/test_swarm.py +++ b/tests/vertexonly/test_swarm.py @@ -26,7 +26,8 @@ def cell_midpoints(m): # MPI rank local. num_cells_local = m.cell_set.size num_cells = MPI.COMM_WORLD.allreduce(num_cells_local, op=MPI.SUM) - local_midpoints = f.dat.data_ro + # reshape is for 1D case where f.dat.data_ro has shape (num_cells_local,) + local_midpoints = f.dat.data_ro.reshape(num_cells_local, m.ufl_cell().geometric_dimension()) local_midpoints_size = np.array(local_midpoints.size) local_midpoints_sizes = np.empty(MPI.COMM_WORLD.size, dtype=int) MPI.COMM_WORLD.Allgatherv(local_midpoints_size, local_midpoints_sizes) @@ -65,9 +66,17 @@ def parentmesh(request): return m +@pytest.fixture(params=["redundant", "nonredundant"]) +def redundant(request): + if request.param == "redundant": + return True + else: + return False + + # pic swarm tests -def test_pic_swarm_in_mesh(parentmesh): +def test_pic_swarm_in_mesh(parentmesh, redundant): """Generate points in cell midpoints of mesh `parentmesh` and check correct swarm is created in plex.""" @@ -78,7 +87,21 @@ def test_pic_swarm_in_mesh(parentmesh): plex = parentmesh.topology.topology_dm from firedrake.petsc import PETSc fields = [("fieldA", 1, PETSc.IntType), ("fieldB", 2, PETSc.ScalarType)] - swarm = mesh._pic_swarm_in_mesh(parentmesh, inputpointcoords, fields=fields) + + if redundant: + # check redundant argument broadcasts from rank 0 by only supplying the + # global cell midpoints only on rank 0. Note that this is the default + # behaviour so it needn't be specified explicitly. + if MPI.COMM_WORLD.rank == 0: + swarm = mesh._pic_swarm_in_mesh(parentmesh, inputpointcoords, fields=fields) + else: + swarm = mesh._pic_swarm_in_mesh(parentmesh, np.empty(inputpointcoords.shape), fields=fields) + else: + # When redundant == False we expect the same behaviour by only + # supplying the local cell midpoints on each MPI ranks. Note that this + # is not the default behaviour so it must be specified explicitly. + swarm = mesh._pic_swarm_in_mesh(parentmesh, inputlocalpointcoords, fields=fields, redundant=redundant) + # Get point coords on current MPI rank localpointcoords = np.copy(swarm.getField("DMSwarmPIC_coor")) swarm.restoreField("DMSwarmPIC_coor") @@ -138,15 +161,17 @@ def test_pic_swarm_in_mesh(parentmesh): @pytest.mark.parallel -def test_pic_swarm_in_mesh_parallel(parentmesh): - test_pic_swarm_in_mesh(parentmesh) +def test_pic_swarm_in_mesh_parallel(parentmesh, redundant): + test_pic_swarm_in_mesh(parentmesh, redundant) @pytest.mark.parallel(nprocs=2) # nprocs == total number of mesh cells def test_pic_swarm_in_mesh_2d_2procs(): - test_pic_swarm_in_mesh(UnitSquareMesh(1, 1)) + test_pic_swarm_in_mesh(UnitSquareMesh(1, 1), redundant=False) + test_pic_swarm_in_mesh(UnitSquareMesh(1, 1), redundant=True) @pytest.mark.parallel(nprocs=3) # nprocs > total number of mesh cells def test_pic_swarm_in_mesh_2d_3procs(): - test_pic_swarm_in_mesh(UnitSquareMesh(1, 1)) + test_pic_swarm_in_mesh(UnitSquareMesh(1, 1), redundant=False) + test_pic_swarm_in_mesh(UnitSquareMesh(1, 1), redundant=False) diff --git a/tests/vertexonly/test_vertex_only_mesh_generation.py b/tests/vertexonly/test_vertex_only_mesh_generation.py index aed69dfb51..14f440386a 100644 --- a/tests/vertexonly/test_vertex_only_mesh_generation.py +++ b/tests/vertexonly/test_vertex_only_mesh_generation.py @@ -26,7 +26,8 @@ def cell_midpoints(m): # MPI rank local. num_cells_local = m.cell_set.size num_cells = MPI.COMM_WORLD.allreduce(num_cells_local, op=MPI.SUM) - local_midpoints = f.dat.data_ro + # reshape is for 1D case where f.dat.data_ro has shape (num_cells_local,) + local_midpoints = f.dat.data_ro.reshape(num_cells_local, m.ufl_cell().geometric_dimension()) local_midpoints_size = np.array(local_midpoints.size) local_midpoints_sizes = np.empty(MPI.COMM_WORLD.size, dtype=int) MPI.COMM_WORLD.Allgatherv(local_midpoints_size, local_midpoints_sizes) @@ -65,6 +66,14 @@ def parentmesh(request): return m +@pytest.fixture(params=["redundant", "nonredundant"]) +def redundant(request): + if request.param == "redundant": + return True + else: + return False + + @pytest.fixture(params=[0, 1, 100], ids=lambda x: f"{x}-coords") def vertexcoords(request, parentmesh): size = (request.param, parentmesh.geometric_dimension()) @@ -133,12 +142,26 @@ def verify_vertexonly_mesh(m, vm, inputvertexcoords): assert m.locate_cell(stored_vertex_coords[i]) == stored_parent_cell_nums[i] -def test_generate_cell_midpoints(parentmesh): +def test_generate_cell_midpoints(parentmesh, redundant): """ Generate cell midpoints for mesh parentmesh and check they lie in the correct cells """ inputcoords, inputcoordslocal = cell_midpoints(parentmesh) + if redundant: + # check redundant argument broadcasts from rank 0 by only supplying the + # global cell midpoints only on rank 0. Note that this is the default + # behaviour so it needn't be specified explicitly. + if MPI.COMM_WORLD.rank == 0: + vm = VertexOnlyMesh(parentmesh, inputcoords) + else: + vm = VertexOnlyMesh(parentmesh, np.empty(inputcoords.shape)) + else: + # When redundant == False we expect the same behaviour by only + # supplying the local cell midpoints on each MPI ranks. Note that this + # is not the default behaviour so it must be specified explicitly. + vm = VertexOnlyMesh(parentmesh, inputcoordslocal, redundant=False) + vm = VertexOnlyMesh(parentmesh, inputcoords) # Midpoints located in correct cells of parent mesh V = VectorFunctionSpace(parentmesh, "DG", 0) @@ -157,8 +180,8 @@ def test_generate_cell_midpoints(parentmesh): @pytest.mark.parallel -def test_generate_cell_midpoints_parallel(parentmesh): - test_generate_cell_midpoints(parentmesh) +def test_generate_cell_midpoints_parallel(parentmesh, redundant): + test_generate_cell_midpoints(parentmesh, redundant) def test_generate_random(parentmesh, vertexcoords): From c0b45ce2123fdeadf358df1d5655ce42f3b3d74b Mon Sep 17 00:00:00 2001 From: Uygar <83709209+uygarkov@users.noreply.github.com> Date: Wed, 1 Feb 2023 17:16:40 +0100 Subject: [PATCH 13/14] Rename .split() -> .subfunctions (#2725) .split() is now marked as deprecated and will spark a warning when used (but will still work) --- demos/camassa-holm/camassaholm.py.rst | 10 +++--- .../linear_fluid_structure_interaction.py.rst | 2 +- demos/ma-demo/ma-demo.py.rst | 2 +- demos/matrix_free/navier_stokes.py.rst | 2 +- demos/matrix_free/rayleigh-benard.py.rst | 2 +- demos/matrix_free/stokes.py.rst | 2 +- demos/multigrid/geometric_multigrid.py.rst | 2 +- demos/poisson/poisson_mixed.py.rst | 2 +- docs/notebooks/05-mixed-poisson.ipynb | 2 +- .../06-pde-constrained-optimisation.ipynb | 4 +-- docs/notebooks/07-geometric-multigrid.ipynb | 2 +- docs/notebooks/08-composable-solvers.ipynb | 2 +- firedrake/adjoint/blocks.py | 2 +- firedrake/assemble.py | 2 +- firedrake/bcs.py | 2 +- firedrake/checkpointing.py | 2 +- firedrake/constant.py | 8 ++++- firedrake/ensemble.py | 14 ++++---- firedrake/formmanipulation.py | 2 +- firedrake/function.py | 36 ++++++++++--------- firedrake/functionspaceimpl.py | 32 ++++++++++++----- firedrake/mg/embedded.py | 4 +-- firedrake/mg/interface.py | 6 ++-- firedrake/preconditioners/patch.py | 6 ++-- firedrake/preconditioners/pmg.py | 2 +- firedrake/slate/slac/kernel_builder.py | 4 +-- firedrake/slate/slate.py | 6 ++-- .../static_condensation/hybridization.py | 20 +++++------ .../slate/static_condensation/la_utils.py | 4 +-- firedrake/slate/static_condensation/scpc.py | 8 ++--- firedrake/solving_utils.py | 2 +- firedrake/tsfc_interface.py | 2 +- firedrake/ufl_expr.py | 6 ++-- tests/equation_bcs/test_bcs_reconstruct.py | 4 +-- tests/equation_bcs/test_equation_bcs.py | 2 +- tests/extrusion/test_mixed_bcs.py | 6 ++-- tests/extrusion/test_real_tensorproduct.py | 2 +- tests/multigrid/test_nested_split.py | 4 +-- tests/multigrid/test_poisson_gtmg.py | 4 +-- tests/multigrid/test_two_poisson_gmg.py | 4 +-- .../output/test_adjoint_disk_checkpointing.py | 2 +- tests/output/test_io_function.py | 4 +-- tests/regression/test_assemble.py | 6 ++-- tests/regression/test_auxiliary_dm.py | 2 +- tests/regression/test_ensembleparallelism.py | 2 +- tests/regression/test_expressions.py | 6 ++-- .../regression/test_fieldsplit_breadcrumbs.py | 2 +- ...est_fieldsplit_fieldsplit_aux_multigrid.py | 10 +++--- .../test_fieldsplit_split_reorder_bcs.py | 2 +- tests/regression/test_function_spaces.py | 4 +-- tests/regression/test_helmholtz_sphere.py | 2 +- tests/regression/test_l2pullbacks.py | 2 +- tests/regression/test_manifolds.py | 4 +-- .../regression/test_mixed_interior_facets.py | 2 +- tests/regression/test_mixed_tensor.py | 2 +- tests/regression/test_moore_spence.py | 6 ++-- tests/regression/test_mtw.py | 2 +- .../regression/test_nonlinear_stokes_hdiv.py | 2 +- tests/regression/test_nullspace.py | 10 +++--- tests/regression/test_par_loops.py | 4 +-- tests/regression/test_patch_pc.py | 2 +- tests/regression/test_piola_mixed_fn.py | 2 +- tests/regression/test_point_eval_api.py | 2 +- tests/regression/test_point_eval_fs.py | 2 +- tests/regression/test_poisson_mixed_no_bcs.py | 2 +- .../test_poisson_mixed_strong_bcs.py | 2 +- tests/regression/test_poisson_sphere.py | 2 +- tests/regression/test_projection.py | 2 +- tests/regression/test_real_space.py | 6 ++-- tests/regression/test_split.py | 4 +-- tests/regression/test_stokes_hdiv_parallel.py | 2 +- tests/regression/test_stokes_mini.py | 2 +- .../test_vector_laplace_on_quadrilaterals.py | 2 +- tests/regression/test_vfs_component_bcs.py | 4 +-- tests/slate/test_assemble_tensors.py | 6 ++-- tests/slate/test_darcy_hybridized_mixed.py | 4 +-- tests/slate/test_hdg_poisson.py | 2 +- tests/slate/test_hybrid_poisson_sphere.py | 2 +- tests/slate/test_matrix_free_hybridization.py | 4 +-- tests/slate/test_slate_hybridization.py | 26 +++++++------- tests/slate/test_slate_hybridization_extr.py | 4 +-- .../slate/test_slate_hybridized_mixed_bcs.py | 8 ++--- tests/slate/test_slate_mixed_direct.py | 2 +- tests/slate/test_variational_prb.py | 2 +- .../test_interpolation_from_parent.py | 6 ++-- 85 files changed, 216 insertions(+), 192 deletions(-) diff --git a/demos/camassa-holm/camassaholm.py.rst b/demos/camassa-holm/camassaholm.py.rst index 6d4330ea2b..740ef24c5e 100644 --- a/demos/camassa-holm/camassaholm.py.rst +++ b/demos/camassa-holm/camassaholm.py.rst @@ -94,11 +94,11 @@ two variables. :: W = MixedFunctionSpace((V, V)) We construct a :class:`~.Function` to store the two variables at time -level ``n``, and :meth:`~.Function.split` it so that we can +level ``n``, and :attr:`~.Function.subfunctions` it so that we can interpolate the initial condition into the two components. :: w0 = Function(W) - m0, u0 = w0.split() + m0, u0 = w0.subfunctions Then we interpolate the initial condition, @@ -161,12 +161,12 @@ rather than blocked system. :: 'ksp_type': 'preonly', 'pc_type': 'lu'}) -Next we use the other form of :meth:`~.Function.split`, ``w0.split()``, +Next we use the other form of :attr:`~.Function.subfunctions`, ``w0.subfunctions``, which is the way to split up a Function in order to access its data e.g. for output. :: - m0, u0 = w0.split() - m1, u1 = w1.split() + m0, u0 = w0.subfunctions + m1, u1 = w1.subfunctions We choose a final time, and initialise a :class:`~.File` object for storing ``u``. as well as an array for storing the function to be visualised:: diff --git a/demos/linear_fluid_structure_interaction/linear_fluid_structure_interaction.py.rst b/demos/linear_fluid_structure_interaction/linear_fluid_structure_interaction.py.rst index 8615755495..015ae52f0f 100644 --- a/demos/linear_fluid_structure_interaction/linear_fluid_structure_interaction.py.rst +++ b/demos/linear_fluid_structure_interaction/linear_fluid_structure_interaction.py.rst @@ -306,7 +306,7 @@ In the end, we proceed with the actual computation loop:: # symplectic Euler scheme LVS_phi_f.solve() LVS_U_phi.solve() - tmp_f, tmp_s = result_mixed.split() + tmp_f, tmp_s = result_mixed.subfunctions phi.assign(tmp_f) U.assign(tmp_s) LVS_eta.solve() diff --git a/demos/ma-demo/ma-demo.py.rst b/demos/ma-demo/ma-demo.py.rst index 9b9a64e3e7..f4b093d5b3 100644 --- a/demos/ma-demo/ma-demo.py.rst +++ b/demos/ma-demo/ma-demo.py.rst @@ -186,7 +186,7 @@ We then put all of these options into the iterative solver, :: and output the solution to a file. :: - u, sigma = w.split() + u, sigma = w.subfunctions u_solv.solve() File("u.pvd").write(u) diff --git a/demos/matrix_free/navier_stokes.py.rst b/demos/matrix_free/navier_stokes.py.rst index eb89b746ab..c74afe3d32 100644 --- a/demos/matrix_free/navier_stokes.py.rst +++ b/demos/matrix_free/navier_stokes.py.rst @@ -130,7 +130,7 @@ find the Reynolds number. :: And finally we write the results to a file for visualisation. :: - u, p = up.split() + u, p = up.subfunctions u.rename("Velocity") p.rename("Pressure") diff --git a/demos/matrix_free/rayleigh-benard.py.rst b/demos/matrix_free/rayleigh-benard.py.rst index 6034371df4..572b5f861d 100644 --- a/demos/matrix_free/rayleigh-benard.py.rst +++ b/demos/matrix_free/rayleigh-benard.py.rst @@ -236,7 +236,7 @@ them, although doing so would be quite easy.:: Finally, we'll output the results for visualisation. :: - u, p, T = upT.split() + u, p, T = upT.subfunctions u.rename("Velocity") p.rename("Pressure") T.rename("Temperature") diff --git a/demos/matrix_free/stokes.py.rst b/demos/matrix_free/stokes.py.rst index 7d186f8d7d..5817ddca6b 100644 --- a/demos/matrix_free/stokes.py.rst +++ b/demos/matrix_free/stokes.py.rst @@ -126,7 +126,7 @@ visualisation. We split the function into its velocity and pressure parts and give them reasonable names, then write them to a paraview file.:: - u, p = up.split() + u, p = up.subfunctions u.rename("Velocity") p.rename("Pressure") diff --git a/demos/multigrid/geometric_multigrid.py.rst b/demos/multigrid/geometric_multigrid.py.rst index abc4ff2118..95759c3184 100644 --- a/demos/multigrid/geometric_multigrid.py.rst +++ b/demos/multigrid/geometric_multigrid.py.rst @@ -255,7 +255,7 @@ approximations. Finally, we'll write the solution for visualisation with Paraview. :: - u, p = u.split() + u, p = u.subfunctions u.rename("Velocity") p.rename("Pressure") diff --git a/demos/poisson/poisson_mixed.py.rst b/demos/poisson/poisson_mixed.py.rst index a3d9b7e7b7..1fafb62d4c 100644 --- a/demos/poisson/poisson_mixed.py.rst +++ b/demos/poisson/poisson_mixed.py.rst @@ -138,7 +138,7 @@ solver parameters. Afterwards we extract the components ``sigma`` and ``u`` on each of the subspaces with ``split``. :: solve(a == L, w, bcs=[bc0, bc1]) - sigma, u = w.split() + sigma, u = w.subfunctions Lastly we write the component of the solution corresponding to the primal variable on the DG space to a file in VTK format for later inspection with a diff --git a/docs/notebooks/05-mixed-poisson.ipynb b/docs/notebooks/05-mixed-poisson.ipynb index 16a43c9e38..d73b536ad9 100644 --- a/docs/notebooks/05-mixed-poisson.ipynb +++ b/docs/notebooks/05-mixed-poisson.ipynb @@ -1755,7 +1755,7 @@ ], "source": [ "# NBVAL_IGNORE_OUTPUT\n", - "sigmah, uh = wh.split()\n", + "sigmah, uh = wh.subfunctions\n", "fig, axes = plt.subplots(ncols=2, sharex=True, sharey=True)\n", "\n", "quiver(sigmah, axes=axes[0])\n", diff --git a/docs/notebooks/06-pde-constrained-optimisation.ipynb b/docs/notebooks/06-pde-constrained-optimisation.ipynb index 6520b8e7ce..824e4c6119 100644 --- a/docs/notebooks/06-pde-constrained-optimisation.ipynb +++ b/docs/notebooks/06-pde-constrained-optimisation.ipynb @@ -2194,7 +2194,7 @@ ], "source": [ "# NBVAL_IGNORE_OUTPUT\n", - "u_init, p_init = w.split()\n", + "u_init, p_init = w.subfunctions\n", "\n", "fig, axes = plt.subplots(nrows=2, sharex=True, sharey=True)\n", "streamlines = streamplot(u_init, resolution=1/3, seed=0, axes=axes[0])\n", @@ -4330,7 +4330,7 @@ ], "source": [ "# NBVAL_IGNORE_OUTPUT\n", - "u_opt, p_opt = w_opt.split()\n", + "u_opt, p_opt = w_opt.subfunctions\n", "\n", "fig, axes = plt.subplots(nrows=2, sharex=True, sharey=True)\n", "streamlines = streamplot(u_opt, resolution=1/3, seed=0, axes=axes[0])\n", diff --git a/docs/notebooks/07-geometric-multigrid.ipynb b/docs/notebooks/07-geometric-multigrid.ipynb index 088a554725..3c4ad69a40 100644 --- a/docs/notebooks/07-geometric-multigrid.ipynb +++ b/docs/notebooks/07-geometric-multigrid.ipynb @@ -1045,7 +1045,7 @@ "source": [ "# NBVAL_IGNORE_OUTPUT\n", "w = solver._problem.u\n", - "u, p = w.split()\n", + "u, p = w.subfunctions\n", "fig, axes = plt.subplots(ncols=2, sharex=True, sharey=True)\n", "streamlines = streamplot(u, resolution=1/30, seed=4, axes=axes[0])\n", "axes[0].set_aspect(\"equal\")\n", diff --git a/docs/notebooks/08-composable-solvers.ipynb b/docs/notebooks/08-composable-solvers.ipynb index ba8f4fb2ea..a0ecbbd89c 100644 --- a/docs/notebooks/08-composable-solvers.ipynb +++ b/docs/notebooks/08-composable-solvers.ipynb @@ -1018,7 +1018,7 @@ "source": [ "%matplotlib notebook\n", "import matplotlib.pyplot as plt\n", - "u_h, p_h = w.split()\n", + "u_h, p_h = w.subfunctions\n", "fig, axes = plt.subplots()\n", "streamlines = streamplot(u_h, resolution=1/30, seed=0, axes=axes)\n", "fig.colorbar(streamlines);" diff --git a/firedrake/adjoint/blocks.py b/firedrake/adjoint/blocks.py index 6018ead99a..967b33bbe5 100644 --- a/firedrake/adjoint/blocks.py +++ b/firedrake/adjoint/blocks.py @@ -321,7 +321,7 @@ def __init__(self, func, idx, ad_block_tag=None): def evaluate_adj_component(self, inputs, adj_inputs, block_variable, idx, prepared=None): if idx == 0: - return adj_inputs[0].split()[self.idx].vector() + return adj_inputs[0].subfunctions[self.idx].vector() else: return adj_inputs[0] diff --git a/firedrake/assemble.py b/firedrake/assemble.py index de1aec7792..9af37a0898 100644 --- a/firedrake/assemble.py +++ b/firedrake/assemble.py @@ -1211,7 +1211,7 @@ def iter_active_coefficients(form, kinfo): """Yield the form coefficients referenced in ``kinfo``.""" for idx, subidxs in kinfo.coefficient_map: for subidx in subidxs: - yield form.coefficients()[idx].split()[subidx] + yield form.coefficients()[idx].subfunctions[subidx] @staticmethod def index_function_spaces(form, indices): diff --git a/firedrake/bcs.py b/firedrake/bcs.py index c660585c5a..7dea14c765 100644 --- a/firedrake/bcs.py +++ b/firedrake/bcs.py @@ -218,7 +218,7 @@ def as_subspace(self, field, V, use_split): if len(field) == 1: W = V else: - W = V.split()[field_renumbering[index]] if use_split else V.sub(field_renumbering[index]) + W = V.subfunctions[field_renumbering[index]] if use_split else V.sub(field_renumbering[index]) if cmpt is not None: W = W.sub(cmpt) return W diff --git a/firedrake/checkpointing.py b/firedrake/checkpointing.py index fd61b1d6ed..ae16126d3b 100644 --- a/firedrake/checkpointing.py +++ b/firedrake/checkpointing.py @@ -780,7 +780,7 @@ def save_function(self, f, idx=None, name=None): if isinstance(V.topological, impl.MixedFunctionSpace): base_path = self._path_to_mixed_function(mesh.name, V_name, f.name()) self.require_group(base_path) - for i, fsub in enumerate(f.split()): + for i, fsub in enumerate(f.subfunctions): path = os.path.join(base_path, str(i)) self.require_group(path) self.set_attr(path, PREFIX + "_function", fsub.name()) diff --git a/firedrake/constant.py b/firedrake/constant.py index fd4d7555f8..fd4342ade9 100644 --- a/firedrake/constant.py +++ b/firedrake/constant.py @@ -111,9 +111,15 @@ def function_space(self): """Return a null function space.""" return None - def split(self): + @utils.cached_property + def subfunctions(self): return (self,) + def split(self): + import warnings + warnings.warn("The .split() method is deprecated, please use the .subfunctions property instead", category=FutureWarning) + return self.subfunctions + def cell_node_map(self, bcs=None): """Return a null cell to node map.""" if bcs is not None: diff --git a/firedrake/ensemble.py b/firedrake/ensemble.py index 61d80a77b9..1496139715 100644 --- a/firedrake/ensemble.py +++ b/firedrake/ensemble.py @@ -100,7 +100,7 @@ def iallreduce(self, f, f_reduced, op=MPI.SUM): :arg f: The a :class:`.Function` to allreduce. :arg f_reduced: the result of the reduction. :arg op: MPI reduction operator. Defaults to MPI.SUM. - :returns: list of MPI.Request objects (one for each of f.split()). + :returns: list of MPI.Request objects (one for each of f.subfunctions). :raises ValueError: if function communicators mismatch each other or the ensemble spatial communicator, or if the functions are in different spaces """ @@ -141,7 +141,7 @@ def ireduce(self, f, f_reduced, op=MPI.SUM, root=0): :arg f_reduced: the result of the reduction on rank root. :arg op: MPI reduction operator. Defaults to MPI.SUM. :arg root: rank to reduce to. Defaults to 0. - :returns: list of MPI.Request objects (one for each of f.split()). + :returns: list of MPI.Request objects (one for each of f.subfunctions). :raises ValueError: if function communicators mismatch each other or the ensemble spatial communicator, or is the functions are in different spaces """ @@ -172,7 +172,7 @@ def ibcast(self, f, root=0): :arg f: The :class:`.Function` to broadcast. :arg root: rank to broadcast from. Defaults to 0. - :returns: list of MPI.Request objects (one for each of f.split()). + :returns: list of MPI.Request objects (one for each of f.subfunctions). :raises ValueError: if function communicator mismatches the ensemble spatial communicator. """ self._check_function(f) @@ -204,7 +204,7 @@ def recv(self, f, source=MPI.ANY_SOURCE, tag=MPI.ANY_TAG, statuses=None): :arg f: The a :class:`.Function` to receive into :arg source: the rank to receive from. Defaults to MPI.ANY_SOURCE. :arg tag: the tag of the message. Defaults to MPI.ANY_TAG. - :arg statuses: MPI.Status objects (one for each of f.split() or None). + :arg statuses: MPI.Status objects (one for each of f.subfunctions or None). :raises ValueError: if function communicator mismatches the ensemble spatial communicator. """ self._check_function(f) @@ -222,7 +222,7 @@ def isend(self, f, dest, tag=0): :arg f: The a :class:`.Function` to send :arg dest: the rank to send to :arg tag: the tag of the message. Defaults to 0. - :returns: list of MPI.Request objects (one for each of f.split()). + :returns: list of MPI.Request objects (one for each of f.subfunctions). :raises ValueError: if function communicator mismatches the ensemble spatial communicator. """ self._check_function(f) @@ -238,7 +238,7 @@ def irecv(self, f, source=MPI.ANY_SOURCE, tag=MPI.ANY_TAG): :arg f: The a :class:`.Function` to receive into :arg source: the rank to receive from. Defaults to MPI.ANY_SOURCE. :arg tag: the tag of the message. Defaults to MPI.ANY_TAG. - :returns: list of MPI.Request objects (one for each of f.split()). + :returns: list of MPI.Request objects (one for each of f.subfunctions). :raises ValueError: if function communicator mismatches the ensemble spatial communicator. """ self._check_function(f) @@ -280,7 +280,7 @@ def isendrecv(self, fsend, dest, sendtag=0, frecv=None, source=MPI.ANY_SOURCE, r :arg frecv: The a :class:`.Function` to receive into. :arg source: the rank to receive from. Defaults to MPI.ANY_SOURCE. :arg recvtag: the tag of the received message. Defaults to MPI.ANY_TAG. - :returns: list of MPI.Request objects (one for each of fsend.split() and frecv.split()). + :returns: list of MPI.Request objects (one for each of fsend.subfunctions and frecv.subfunctions). :raises ValueError: if function communicator mismatches the ensemble spatial communicator. """ # functions don't necessarily have to match diff --git a/firedrake/formmanipulation.py b/firedrake/formmanipulation.py index 431137985c..7b14f6b8a3 100644 --- a/firedrake/formmanipulation.py +++ b/firedrake/formmanipulation.py @@ -95,7 +95,7 @@ def argument(self, o): if o in self._arg_cache: return self._arg_cache[o] - V_is = V.split() + V_is = V.subfunctions indices = self.blocks[o.number()] try: diff --git a/firedrake/function.py b/firedrake/function.py index 1c697763bc..dbc7b7c6ae 100644 --- a/firedrake/function.py +++ b/firedrake/function.py @@ -106,16 +106,18 @@ def ufl_id(self): return self.uid @utils.cached_property - def _split(self): + def subfunctions(self): + r"""Extract any sub :class:`Function`\s defined on the component spaces + of this this :class:`Function`'s :class:`.FunctionSpace`.""" return tuple(CoordinatelessFunction(fs, dat, name="%s[%d]" % (self.name(), i)) for i, (fs, dat) in enumerate(zip(self.function_space(), self.dat))) @PETSc.Log.EventDecorator() def split(self): - r"""Extract any sub :class:`Function`\s defined on the component spaces - of this this :class:`Function`'s :class:`.FunctionSpace`.""" - return self._split + import warnings + warnings.warn("The .split() method is deprecated, please use the .subfunctions property instead", category=FutureWarning) + return self.subfunctions @utils.cached_property def _components(self): @@ -132,7 +134,7 @@ def sub(self, i): :arg i: the index to extract - See also :meth:`split`. + See also :attr:`subfunctions`. If the :class:`Function` is defined on a rank-n :class:`~.FunctionSpace`, this returns a proxy object @@ -140,7 +142,7 @@ def sub(self, i): boundary condition application.""" if len(self.function_space()) == 1: return self._components[i] - return self._split[i] + return self.subfunctions[i] @property def cell_set(self): @@ -302,16 +304,18 @@ def __dir__(self): return list(OrderedDict.fromkeys(dir(self._data) + current)) @utils.cached_property - def _split(self): + def subfunctions(self): + r"""Extract any sub :class:`Function`\s defined on the component spaces + of this this :class:`Function`'s :class:`.FunctionSpace`.""" return tuple(type(self)(V, val) - for (V, val) in zip(self.function_space(), self.topological.split())) + for (V, val) in zip(self.function_space(), self.topological.subfunctions)) @PETSc.Log.EventDecorator() @FunctionMixin._ad_annotate_split def split(self): - r"""Extract any sub :class:`Function`\s defined on the component spaces - of this this :class:`Function`'s :class:`.FunctionSpace`.""" - return self._split + import warnings + warnings.warn("The .split() method is deprecated, please use the .subfunctions property instead", category=FutureWarning) + return self.subfunctions @utils.cached_property def _components(self): @@ -327,7 +331,7 @@ def sub(self, i): :arg i: the index to extract - See also :meth:`split`. + See also :attr:`subfunctions`. If the :class:`Function` is defined on a :class:`~.VectorFunctionSpace` or :class:`~.TensorFunctiionSpace` this returns a proxy object @@ -335,7 +339,7 @@ def sub(self, i): boundary condition application.""" if len(self.function_space()) == 1: return self._components[i] - return self._split[i] + return self.subfunctions[i] @PETSc.Log.EventDecorator() @FunctionMixin._ad_annotate_project @@ -549,15 +553,15 @@ def single_eval(x, buf): points = arg.reshape(-1, arg.shape[-1]) value_shape = self.ufl_shape - split = self.split() - mixed = len(split) != 1 + subfunctions = self.subfunctions + mixed = len(subfunctions) != 1 # Local evaluation l_result = [] for i, p in enumerate(points): try: if mixed: - l_result.append((i, tuple(f.at(p) for f in split))) + l_result.append((i, tuple(f.at(p) for f in subfunctions))) else: p_result = np.zeros(value_shape, dtype=ScalarType) single_eval(points[i:i+1], p_result) diff --git a/firedrake/functionspaceimpl.py b/firedrake/functionspaceimpl.py index f91ed54983..048391d555 100644 --- a/firedrake/functionspaceimpl.py +++ b/firedrake/functionspaceimpl.py @@ -100,9 +100,10 @@ def topological(self, val): self.cargo.topological = val @utils.cached_property - def _split(self): + def subfunctions(self): + r"""Split into a tuple of constituent spaces.""" return tuple(WithGeometry.create(subspace, self.mesh()) - for subspace in self.topological.split()) + for subspace in self.topological.subfunctions) mesh = ufl.FunctionSpace.ufl_domain @@ -120,8 +121,9 @@ def ufl_cell(self): @PETSc.Log.EventDecorator() def split(self): - r"""Split into a tuple of constituent spaces.""" - return self._split + import warnings + warnings.warn("The .split() method is deprecated, please use the .subfunctions property instead", category=FutureWarning) + return self.subfunctions @utils.cached_property def _components(self): @@ -129,7 +131,7 @@ def _components(self): return tuple(WithGeometry.create(self.topological.sub(i), self.mesh()) for i in range(self.value_size)) else: - return self._split + return self.subfunctions @PETSc.Log.EventDecorator() def sub(self, i): @@ -269,10 +271,10 @@ def __str__(self): return "WithGeometry(%s, %s)" % (self.topological, self.mesh()) def __iter__(self): - return iter(self._split) + return iter(self.subfunctions) def __getitem__(self, i): - return self._split[i] + return self.subfunctions[i] def __mul__(self, other): r"""Create a :class:`.MixedFunctionSpace` composed of this @@ -496,10 +498,16 @@ def __str__(self): self.ufl_element(), self.name) - def split(self): + @utils.cached_property + def subfunctions(self): r"""Split into a tuple of constituent spaces.""" return (self, ) + def split(self): + import warnings + warnings.warn("The .split() method is deprecated, please use the .subfunctions property instead", category=FutureWarning) + return self.subfunctions + def __getitem__(self, i): r"""Return the ith subspace.""" if i != 0: @@ -705,11 +713,17 @@ def __ne__(self, other): def __hash__(self): return hash(tuple(self)) - def split(self): + @utils.cached_property + def subfunctions(self): r"""The list of :class:`FunctionSpace`\s of which this :class:`MixedFunctionSpace` is composed.""" return self._spaces + def split(self): + import warnings + warnings.warn("The .split() method is deprecated, please use the .subfunctions property instead", category=FutureWarning) + return self.subfunctions + def sub(self, i): r"""Return the `i`th :class:`FunctionSpace` in this :class:`MixedFunctionSpace`.""" diff --git a/firedrake/mg/embedded.py b/firedrake/mg/embedded.py index c7d977818b..e719c3615d 100644 --- a/firedrake/mg/embedded.py +++ b/firedrake/mg/embedded.py @@ -211,7 +211,7 @@ def op(self, source, target, transfer_op): return self._native_transfer(source_element, transfer_op)(source, target) if type(source_element) is ufl.MixedElement: assert type(target_element) is ufl.MixedElement - for source_, target_ in zip(source.split(), target.split()): + for source_, target_ in zip(source.subfunctions, target.subfunctions): self.op(source_, target_, transfer_op=transfer_op) return target # Get some work vectors @@ -273,7 +273,7 @@ def restrict(self, gf, gc): return self._native_transfer(source_element, Op.RESTRICT)(gf, gc) if type(source_element) is ufl.MixedElement: assert type(target_element) is ufl.MixedElement - for source_, target_ in zip(gf.split(), gc.split()): + for source_, target_ in zip(gf.subfunctions, gc.subfunctions): self.restrict(source_, target_) return gc dgf = self.DG_work(Vf) diff --git a/firedrake/mg/interface.py b/firedrake/mg/interface.py index 0ce9d2af00..8d6d5bc7c7 100644 --- a/firedrake/mg/interface.py +++ b/firedrake/mg/interface.py @@ -32,7 +32,7 @@ def prolong(coarse, fine): if len(Vc) > 1: if len(Vc) != len(Vf): raise ValueError("Mixed spaces have different lengths") - for in_, out in zip(coarse.split(), fine.split()): + for in_, out in zip(coarse.subfunctions, fine.subfunctions): manager = firedrake.dmhooks.get_transfer_manager(in_.function_space().dm) manager.prolong(in_, out) return fine @@ -93,7 +93,7 @@ def restrict(fine_dual, coarse_dual): if len(Vc) > 1: if len(Vc) != len(Vf): raise ValueError("Mixed spaces have different lengths") - for in_, out in zip(fine_dual.split(), coarse_dual.split()): + for in_, out in zip(fine_dual.subfunctions, coarse_dual.subfunctions): manager = firedrake.dmhooks.get_transfer_manager(in_.function_space().dm) manager.restrict(in_, out) return coarse_dual @@ -155,7 +155,7 @@ def inject(fine, coarse): if len(Vc) > 1: if len(Vc) != len(Vf): raise ValueError("Mixed spaces have different lengths") - for in_, out in zip(fine.split(), coarse.split()): + for in_, out in zip(fine.subfunctions, coarse.subfunctions): manager = firedrake.dmhooks.get_transfer_manager(in_.function_space().dm) manager.inject(in_, out) return diff --git a/firedrake/preconditioners/patch.py b/firedrake/preconditioners/patch.py index e2d346838e..2a20f2e2f8 100644 --- a/firedrake/preconditioners/patch.py +++ b/firedrake/preconditioners/patch.py @@ -205,7 +205,7 @@ def matrix_funptr(form, state): args.append(statearg) continue for ind in indices: - c_ = c.split()[ind] + c_ = c.subfunctions[ind] map_ = get_map(c_) arg = c_.dat(op2.READ, map_) args.append(arg) @@ -295,7 +295,7 @@ def residual_funptr(form, state): args.append(statearg) continue for ind in indices: - c_ = c.split()[ind] + c_ = c.subfunctions[ind] map_ = get_map(c_) arg = c_.dat(op2.READ, map_) args.append(arg) @@ -516,7 +516,7 @@ def make_c_arguments(form, kernel, state, get_map, require_state=False, raise ValueError(f"Active indices of state (dont_split) function must be (0, ), not {indices}") coeffs.append(c) else: - coeffs.extend([c.split()[ind] for ind in indices]) + coeffs.extend([c.subfunctions[ind] for ind in indices]) if require_state: assert state in coeffs, "Couldn't find state vector in form coefficients" data_args = [] diff --git a/firedrake/preconditioners/pmg.py b/firedrake/preconditioners/pmg.py index dddc26068c..88cc0f1b10 100644 --- a/firedrake/preconditioners/pmg.py +++ b/firedrake/preconditioners/pmg.py @@ -1153,7 +1153,7 @@ def __init__(self, Vf, Vc, Vf_bcs, Vc_bcs): self.uc = Vc if isinstance(Vc, firedrake.Function) else firedrake.Function(Vc) self.standalones = [] - for (i, (uf_sub, uc_sub)) in enumerate(zip(self.uf.split(), self.uc.split())): + for (i, (uf_sub, uc_sub)) in enumerate(zip(self.uf.subfunctions, self.uc.subfunctions)): Vf_sub_bcs = [bc for bc in Vf_bcs if bc.function_space().index == i] Vc_sub_bcs = [bc for bc in Vc_bcs if bc.function_space().index == i] standalone = StandaloneInterpolationMatrix(uf_sub, uc_sub, Vf_sub_bcs, Vc_sub_bcs) diff --git a/firedrake/slate/slac/kernel_builder.py b/firedrake/slate/slac/kernel_builder.py index e53c933aed..fe2f1d60aa 100644 --- a/firedrake/slate/slac/kernel_builder.py +++ b/firedrake/slate/slac/kernel_builder.py @@ -334,7 +334,7 @@ def coefficient_map(self): for i, coefficient in enumerate(self.expression.coefficients()): if type(coefficient.ufl_element()) == MixedElement: csym_info = [] - for j, _ in enumerate(coefficient.split()): + for j, _ in enumerate(coefficient.subfunctions): csym_info.append(ast.Symbol("w_%d_%d" % (i, j))) else: csym_info = (ast.Symbol("w_%d" % i),) @@ -589,7 +589,7 @@ def collect_coefficients(self): for i, (c, split_map) in enumerate(self.expression.coeff_map): coeff = coeffs[c] if type(coeff.ufl_element()) == MixedElement: - splits = coeff.split() + splits = coeff.subfunctions coeff_dict[coeff] = OrderedDict({splits[j]: (f"w_{i}_{j}", self.extent(splits[j])) for j in split_map}) else: coeff_dict[coeff] = (f"w_{i}", self.extent(coeff)) diff --git a/firedrake/slate/slate.py b/firedrake/slate/slate.py index c70bc35870..03311c47d8 100644 --- a/firedrake/slate/slate.py +++ b/firedrake/slate/slate.py @@ -230,7 +230,7 @@ def coeff_map(self): coeff_map[m].update(c.indices[0]) else: m = self.coefficients().index(c) - split_map = tuple(range(len(c.split()))) if isinstance(c, Function) or isinstance(c, Constant) else tuple(range(1)) + split_map = tuple(range(len(c.subfunctions))) if isinstance(c, Function) or isinstance(c, Constant) else tuple(range(1)) coeff_map[m].update(split_map) return tuple((k, tuple(sorted(v)))for k, v in coeff_map.items()) @@ -662,7 +662,7 @@ def _split_arguments(self): nargs = [] for i, arg in enumerate(tensor.arguments()): V = arg.function_space() - V_is = V.split() + V_is = V.subfunctions idx = as_tuple(self._blocks[i]) if len(idx) == 1: fidx, = idx @@ -696,7 +696,7 @@ def form(self): else: # turns the Block on an AssembledVector into a set off coefficients # corresponding to the indices of the Block - return tuple(tensor._function.split()[i] for i in chain(*self._indices)) + return tuple(tensor._function.subfunctions[i] for i in chain(*self._indices)) @cached_property def assembled(self): diff --git a/firedrake/slate/static_condensation/hybridization.py b/firedrake/slate/static_condensation/hybridization.py index 992727d4d9..3b97e2cf7c 100644 --- a/firedrake/slate/static_condensation/hybridization.py +++ b/firedrake/slate/static_condensation/hybridization.py @@ -284,8 +284,8 @@ def _reconstruction_calls(self): S = self.schur_builder.inner_S # Split functions and reconstruct each bit separately - split_residual = self.broken_residual.split() - split_sol = self.broken_solution.split() + split_residual = self.broken_residual.subfunctions + split_sol = self.broken_solution.subfunctions g = AssembledVector(split_residual[id0]) f = AssembledVector(split_residual[id1]) sigma = split_sol[id0] @@ -334,8 +334,8 @@ def forward_elimination(self, pc, x): # Transfer unbroken_rhs into broken_rhs # NOTE: Scalar space is already "broken" so no need for # any projections - unbroken_scalar_data = self.unbroken_residual.split()[self.pidx] - broken_scalar_data = self.broken_residual.split()[self.pidx] + unbroken_scalar_data = self.unbroken_residual.subfunctions[self.pidx] + broken_scalar_data = self.broken_residual.subfunctions[self.pidx] unbroken_scalar_data.dat.copy(broken_scalar_data.dat) # Assemble the new "broken" hdiv residual @@ -344,8 +344,8 @@ def forward_elimination(self, pc, x): # We do this by splitting the residual equally between # basis functions that add together to give unbroken # basis functions. - unbroken_res_hdiv = self.unbroken_residual.split()[self.vidx] - broken_res_hdiv = self.broken_residual.split()[self.vidx] + unbroken_res_hdiv = self.unbroken_residual.subfunctions[self.vidx] + broken_res_hdiv = self.broken_residual.subfunctions[self.vidx] broken_res_hdiv.assign(0) par_loop(self.average_kernel, ufl.dx, {"w": (self.weight, READ), @@ -393,13 +393,13 @@ def backward_substitution(self, pc, y): with PETSc.Log.Event("HybridProject"): # Project the broken solution into non-broken spaces - broken_pressure = self.broken_solution.split()[self.pidx] - unbroken_pressure = self.unbroken_solution.split()[self.pidx] + broken_pressure = self.broken_solution.subfunctions[self.pidx] + unbroken_pressure = self.unbroken_solution.subfunctions[self.pidx] broken_pressure.dat.copy(unbroken_pressure.dat) # Compute the hdiv projection of the broken hdiv solution - broken_hdiv = self.broken_solution.split()[self.vidx] - unbroken_hdiv = self.unbroken_solution.split()[self.vidx] + broken_hdiv = self.broken_solution.subfunctions[self.vidx] + unbroken_hdiv = self.unbroken_solution.subfunctions[self.vidx] unbroken_hdiv.assign(0) par_loop(self.average_kernel, ufl.dx, diff --git a/firedrake/slate/static_condensation/la_utils.py b/firedrake/slate/static_condensation/la_utils.py index f8e8dd5869..96dd84e349 100644 --- a/firedrake/slate/static_condensation/la_utils.py +++ b/firedrake/slate/static_condensation/la_utils.py @@ -119,8 +119,8 @@ def backward_solve(A, b, x, schur_builder, reconstruct_fields): nfields = len(all_fields) reconstruct_fields = as_tuple(reconstruct_fields) - _b = b.split() - _x = x.split() + _b = b.subfunctions + _x = x.subfunctions # Ordering matters systems = [] diff --git a/firedrake/slate/static_condensation/scpc.py b/firedrake/slate/static_condensation/scpc.py index 0afbe8225c..2e6596700f 100644 --- a/firedrake/slate/static_condensation/scpc.py +++ b/firedrake/slate/static_condensation/scpc.py @@ -219,7 +219,7 @@ def local_solver_calls(self, A, rhs, x, elim_fields, schur_builder): from firedrake.assemble import OneFormAssembler from firedrake.slate.static_condensation.la_utils import backward_solve - fields = x.split() + fields = x.subfunctions systems = backward_solve(A, rhs, x, schur_builder, reconstruct_fields=elim_fields) local_solvers = [] @@ -260,7 +260,7 @@ def forward_elimination(self, pc, x): x.copy(v) # Disassemble the incoming right-hand side - with self.residual.split()[self.c_field].dat.vec as vc, self.weight.dat.vec_ro as wc: + with self.residual.subfunctions[self.c_field].dat.vec as vc, self.weight.dat.vec_ro as wc: vc.pointwiseMult(vc, wc) # Now assemble residual for the reduced problem @@ -279,9 +279,9 @@ def sc_solve(self, pc): with self.condensed_rhs.dat.vec_ro as rhs: if self.condensed_ksp.getInitialGuessNonzero(): - acc = self.solution.split()[self.c_field].dat.vec + acc = self.solution.subfunctions[self.c_field].dat.vec else: - acc = self.solution.split()[self.c_field].dat.vec_wo + acc = self.solution.subfunctions[self.c_field].dat.vec_wo with acc as sol: self.condensed_ksp.solve(rhs, sol) diff --git a/firedrake/solving_utils.py b/firedrake/solving_utils.py index 670ffe1b88..5055582f20 100644 --- a/firedrake/solving_utils.py +++ b/firedrake/solving_utils.py @@ -327,7 +327,7 @@ def split(self, fields): for field in fields: F = splitter.split(problem.F, argument_indices=(field, )) J = splitter.split(problem.J, argument_indices=(field, field)) - us = problem.u.split() + us = problem.u.subfunctions V = F.arguments()[0].function_space() # Exposition: # We are going to make a new solution Function on the sub diff --git a/firedrake/tsfc_interface.py b/firedrake/tsfc_interface.py index 10ab7a34f9..bc92def201 100644 --- a/firedrake/tsfc_interface.py +++ b/firedrake/tsfc_interface.py @@ -315,7 +315,7 @@ def extract_numbered_coefficients(expr, numbers): coefficients = [] for coeff in (orig_coefficients[i] for i in numbers): if type(coeff.ufl_element()) == ufl.MixedElement: - coefficients.extend(coeff.split()) + coefficients.extend(coeff.subfunctions) else: coefficients.append(coeff) return coefficients diff --git a/firedrake/ufl_expr.py b/firedrake/ufl_expr.py index 8f187e1f39..d2e86a76b1 100644 --- a/firedrake/ufl_expr.py +++ b/firedrake/ufl_expr.py @@ -135,7 +135,7 @@ def derivative(form, u, du=None, coefficient_derivatives=None): provide the derivative of a coefficient function. :raises ValueError: If any of the coefficients in ``form`` were - obtained from ``u.split()``. UFL doesn't notice that these + obtained from ``u.subfunctions``. UFL doesn't notice that these are related to ``u`` and so therefore the derivative is wrong (instead one should have written ``split(u)``). @@ -148,8 +148,8 @@ def derivative(form, u, du=None, coefficient_derivatives=None): # TODO: What about Constant? u_is_x = isinstance(u, ufl.SpatialCoordinate) uc, = (u,) if u_is_x else extract_coefficients(u) - if not u_is_x and len(uc.split()) > 1 and set(extract_coefficients(form)) & set(uc.split()): - raise ValueError("Taking derivative of form wrt u, but form contains coefficients from u.split()." + if not u_is_x and len(uc.subfunctions) > 1 and set(extract_coefficients(form)) & set(uc.subfunctions): + raise ValueError("Taking derivative of form wrt u, but form contains coefficients from u.subfunctions." "\nYou probably meant to write split(u) when defining your form.") mesh = form.ufl_domain() diff --git a/tests/equation_bcs/test_bcs_reconstruct.py b/tests/equation_bcs/test_bcs_reconstruct.py index 914c993aef..31e6bdfbcf 100755 --- a/tests/equation_bcs/test_bcs_reconstruct.py +++ b/tests/equation_bcs/test_bcs_reconstruct.py @@ -49,5 +49,5 @@ def test_bc_on_sub_sub_domain(): f.interpolate(as_vector([cos(2 * pi * x) * cos(2 * pi * y), cos(2 * pi * x) * cos(2 * pi * y)])) - assert sqrt(assemble(dot(uu.split()[0] - f, uu.split()[0] - f) * dx)) < 4.0e-03 - assert sqrt(assemble(dot(uu.split()[1] - f, uu.split()[1] - f) * dx)) < 4.0e-03 + assert sqrt(assemble(dot(uu.subfunctions[0] - f, uu.subfunctions[0] - f) * dx)) < 4.0e-03 + assert sqrt(assemble(dot(uu.subfunctions[1] - f, uu.subfunctions[1] - f) * dx)) < 4.0e-03 diff --git a/tests/equation_bcs/test_equation_bcs.py b/tests/equation_bcs/test_equation_bcs.py index 45f1674fe5..6dbae85534 100644 --- a/tests/equation_bcs/test_equation_bcs.py +++ b/tests/equation_bcs/test_equation_bcs.py @@ -189,7 +189,7 @@ def linear_poisson_mixed(solver_parameters, mesh_num, porder): solve(a == L, w, bcs=[bc2, bc3, bc4], solver_parameters=solver_parameters) - sigma, u = w.split() + sigma, u = w.subfunctions f.interpolate(cos(2 * pi * x + pi / 3) * cos(2 * pi * y)) g = Function(BDM).project(as_vector([-2 * pi * sin(2 * pi * x + pi / 3) * cos(2 * pi * y), -2 * pi * cos(2 * pi * x + pi / 3) * sin(2 * pi * y)])) diff --git a/tests/extrusion/test_mixed_bcs.py b/tests/extrusion/test_mixed_bcs.py index dcc8ec3d8f..ba6edc92f9 100644 --- a/tests/extrusion/test_mixed_bcs.py +++ b/tests/extrusion/test_mixed_bcs.py @@ -41,7 +41,7 @@ def test_multiple_poisson_Pn(quadrilateral, degree): wexact = Function(W) - u, p = wexact.split() + u, p = wexact.subfunctions xs = SpatialCoordinate(mesh) u.interpolate(1 + 9*xs[2]) @@ -88,7 +88,7 @@ def test_multiple_poisson_strong_weak_Pn(quadrilateral, degree): wexact = Function(W) - u, p = wexact.split() + u, p = wexact.subfunctions xs = SpatialCoordinate(mesh) u.interpolate(11 - xs[2]) @@ -134,7 +134,7 @@ def test_stokes_taylor_hood(mat_type): w = Function(W) - u, p = w.split() + u, p = w.subfunctions solve(a == L, w, bcs=bcs, solver_parameters={'pc_type': 'fieldsplit', 'pc_fieldsplit_type': 'schur', diff --git a/tests/extrusion/test_real_tensorproduct.py b/tests/extrusion/test_real_tensorproduct.py index e9522107a0..2c82d56950 100644 --- a/tests/extrusion/test_real_tensorproduct.py +++ b/tests/extrusion/test_real_tensorproduct.py @@ -167,7 +167,7 @@ def test_real_tensorproduct_mixed(V): Q = FunctionSpace(mesh, "P", 2) W = V*Q - for (s_, s) in zip(W.split(), (V, Q)): + for (s_, s) in zip(W.subfunctions, (V, Q)): assert s_.node_set is s.node_set assert s_.dof_dset is s.dof_dset diff --git a/tests/multigrid/test_nested_split.py b/tests/multigrid/test_nested_split.py index adaa50df2d..80494c3863 100644 --- a/tests/multigrid/test_nested_split.py +++ b/tests/multigrid/test_nested_split.py @@ -85,7 +85,7 @@ def test_nested_split_multigrid(parameters): F += inner(s, r)*dx - inner(x, r)*dx expect = Function(W) - u_expect, p_expect, s_expect = expect.split() + u_expect, p_expect, s_expect = expect.subfunctions u_expect.interpolate(x**2 + y**2) p_expect.interpolate(sin(pi*x)*tan(pi*x*0.25)*sin(pi*y)) @@ -98,7 +98,7 @@ def test_nested_split_multigrid(parameters): solver = NonlinearVariationalSolver(problem, options_prefix="", solver_parameters=parameters) solver.solve() - u, p, s = w.split() + u, p, s = w.subfunctions assert norm(assemble(u_expect - u)) < 5e-5 assert norm(assemble(p_expect - p)) < 1e-6 diff --git a/tests/multigrid/test_poisson_gtmg.py b/tests/multigrid/test_poisson_gtmg.py index 708a69a491..5c2be75b95 100644 --- a/tests/multigrid/test_poisson_gtmg.py +++ b/tests/multigrid/test_poisson_gtmg.py @@ -59,7 +59,7 @@ def p1_callback(): 'coarse_space_bcs': get_p1_prb_bcs()} solve(a == L, w, solver_parameters=params, appctx=appctx) - _, uh = w.split() + _, uh = w.subfunctions # Analytical solution f.interpolate(x[0]*(1-x[0])*x[1]*(1-x[1])) @@ -135,7 +135,7 @@ def p1_callback(): bcs = DirichletBC(W.sub(2), Constant(0.0), "on_boundary") solve(a == L, w, bcs=bcs, solver_parameters=params, appctx=appctx) - _, uh, _ = w.split() + _, uh, _ = w.subfunctions # Analytical solution f.interpolate(x[0]*(1-x[0])*x[1]*(1-x[1])) diff --git a/tests/multigrid/test_two_poisson_gmg.py b/tests/multigrid/test_two_poisson_gmg.py index 1fcaf8800b..966ce2e2b9 100644 --- a/tests/multigrid/test_two_poisson_gmg.py +++ b/tests/multigrid/test_two_poisson_gmg.py @@ -80,7 +80,7 @@ def run_two_poisson(typ): # eigenmode. This stresses the preconditioner much more. e.g. 10 # iterations of ilu fails to converge this problem sufficiently. x = SpatialCoordinate(W.mesh()) - for h in f.split(): + for h in f.subfunctions: h.interpolate(-0.5*pi*pi*(4*cos(pi*x[0]) - 5*cos(pi*x[0]*0.5) + 2)*sin(pi*x[1])) problem = NonlinearVariationalProblem(F, u, bcs=bcs) @@ -94,7 +94,7 @@ def run_two_poisson(typ): for exact in [exact_P2, exact_P1]: exact.interpolate(sin(pi*x[0])*tan(pi*x[0]*0.25)*sin(pi*x[1])) - sol_P2, sol_P1 = u.split() + sol_P2, sol_P1 = u.subfunctions return norm(assemble(exact_P2 - sol_P2)), norm(assemble(exact_P1 - sol_P1)) diff --git a/tests/output/test_adjoint_disk_checkpointing.py b/tests/output/test_adjoint_disk_checkpointing.py index 4543555d81..4c166b25cf 100644 --- a/tests/output/test_adjoint_disk_checkpointing.py +++ b/tests/output/test_adjoint_disk_checkpointing.py @@ -32,7 +32,7 @@ def adjoint_example(mesh): # InterpolateBlock m = interpolate(sin(4*pi*x)*cos(4*pi*y), cg_space) - u, v = w.split() + u, v = w.subfunctions # FunctionAssignBlock, FunctionMergeBlock v.assign(m) # FunctionSplitBlock, GenericSolveBlock diff --git a/tests/output/test_io_function.py b/tests/output/test_io_function.py index 2276a47ea7..588ef6b84b 100644 --- a/tests/output/test_io_function.py +++ b/tests/output/test_io_function.py @@ -263,7 +263,7 @@ def test_io_function_mixed_real(cell_family_degree_tuples, tmpdir): VA = functools.reduce(lambda a, b: a * b, VA_list) method = "project" fA = Function(VA, name=func_name) - fA0, fA1 = fA.split() + fA0, fA1 = fA.subfunctions _initialise_function(fA0, _get_expr(VA[0]), method) fA1.dat.data.itemset(3.14) with CheckpointFile(filename, 'w', comm=COMM_WORLD) as afile: @@ -279,7 +279,7 @@ def test_io_function_mixed_real(cell_family_degree_tuples, tmpdir): fB = afile.load_function(meshB, func_name) VB = fB.function_space() fBe = Function(VB) - fBe0, fBe1 = fBe.split() + fBe0, fBe1 = fBe.subfunctions _initialise_function(fBe0, _get_expr(VB[0]), method) fBe1.dat.data.itemset(3.14) assert assemble(inner(fB - fBe, fB - fBe) * dx) < 1.e-16 diff --git a/tests/regression/test_assemble.py b/tests/regression/test_assemble.py index 2a215c4c03..cd2d2db0d2 100644 --- a/tests/regression/test_assemble.py +++ b/tests/regression/test_assemble.py @@ -41,7 +41,7 @@ def fs(request, mesh): @pytest.fixture def f(fs): f = Function(fs, name="f") - f_split = f.split() + f_split = f.subfunctions x = SpatialCoordinate(fs.mesh())[0] # NOTE: interpolation of UFL expressions into mixed @@ -61,7 +61,7 @@ def f(fs): @pytest.fixture def one(fs): one = Function(fs, name="one") - ones = one.split() + ones = one.subfunctions # NOTE: interpolation of UFL expressions into mixed # function spaces is not yet implemented @@ -87,7 +87,7 @@ def M(fs): def test_one_form(M, f): one_form = assemble(action(M, f)) assert isinstance(one_form, Function) - for f in one_form.split(): + for f in one_form.subfunctions: if f.function_space().rank == 2: assert abs(f.dat.data.sum() - 0.5*sum(f.function_space().shape)) < 1.0e-12 else: diff --git a/tests/regression/test_auxiliary_dm.py b/tests/regression/test_auxiliary_dm.py index 3704f7520a..ac963fa69a 100644 --- a/tests/regression/test_auxiliary_dm.py +++ b/tests/regression/test_auxiliary_dm.py @@ -131,7 +131,7 @@ def test_auxiliary_dm(): solver.solve() # Error in L2 norm - (u, v, alpha) = z.split() + (u, v, alpha) = z.subfunctions u_exact = problem.analytical_solution(mesh) error_L2 = errornorm(u_exact, u, 'L2') / errornorm(u_exact, Function(FunctionSpace(mesh, 'CG', 1)), 'L2') assert error_L2 < 0.02 diff --git a/tests/regression/test_ensembleparallelism.py b/tests/regression/test_ensembleparallelism.py index 66928f2f77..9511b1a47e 100644 --- a/tests/regression/test_ensembleparallelism.py +++ b/tests/regression/test_ensembleparallelism.py @@ -32,7 +32,7 @@ def function_profile(x, y, rank, cpt): def unique_function(mesh, rank, W): u = Function(W) x, y = SpatialCoordinate(mesh) - for cpt, v in enumerate(u.split()): + for cpt, v in enumerate(u.subfunctions): v.interpolate(function_profile(x, y, rank, cpt)) return u diff --git a/tests/regression/test_expressions.py b/tests/regression/test_expressions.py index 4f03959a15..f22c48a22d 100644 --- a/tests/regression/test_expressions.py +++ b/tests/regression/test_expressions.py @@ -374,7 +374,7 @@ def test_assign_from_mfs_sub(cg1, vcg1): u = Function(cg1) v = Function(vcg1) - w1, w2 = w.split() + w1, w2 = w.subfunctions w1.assign(4) w2.assign(10) @@ -389,7 +389,7 @@ def test_assign_from_mfs_sub(cg1, vcg1): Q = vcg1*cg1 q = Function(Q) - q1, q2 = q.split() + q1, q2 = q.subfunctions q1.assign(11) q2.assign(12) @@ -453,7 +453,7 @@ def test_assign_mixed_multiple_shaped(): z2.dat[3].data[:] = [[15, 16], [17, 18]] q = assemble(z1 - z2) - for q, p1, p2 in zip(q.split(), z1.split(), z2.split()): + for q, p1, p2 in zip(q.subfunctions, z1.subfunctions, z2.subfunctions): assert np.allclose(q.dat.data_ro, p1.dat.data_ro - p2.dat.data_ro) diff --git a/tests/regression/test_fieldsplit_breadcrumbs.py b/tests/regression/test_fieldsplit_breadcrumbs.py index 48389553e5..1435586d0e 100644 --- a/tests/regression/test_fieldsplit_breadcrumbs.py +++ b/tests/regression/test_fieldsplit_breadcrumbs.py @@ -34,6 +34,6 @@ def test_fieldsplit_breadcrumbs(): u_source.assign(g*i) solver.solve() - u, eta = solution.split() + u, eta = solution.subfunctions assert numpy.allclose(u.dat.data_ro, u_source.dat.data_ro) assert numpy.allclose(eta.dat.data_ro, 0) diff --git a/tests/regression/test_fieldsplit_fieldsplit_aux_multigrid.py b/tests/regression/test_fieldsplit_fieldsplit_aux_multigrid.py index 881f2d5ae3..f6ebb28e01 100644 --- a/tests/regression/test_fieldsplit_fieldsplit_aux_multigrid.py +++ b/tests/regression/test_fieldsplit_fieldsplit_aux_multigrid.py @@ -94,11 +94,11 @@ def alpha(d): "snes_monitor": None} nsp_guess = MixedVectorSpaceBasis(G, [G.sub(0), VectorSpaceBasis(constant=True)]) solve(F_guess == 0, g, bcs=Gbcs, nullspace=nsp_guess, solver_parameters=solver_params_guess) - (u_guess, p_guess) = g.split() - z.split()[0].interpolate(Constant(gamma)) - z.split()[1].assign(u_guess) - z.split()[2].assign(p_guess) - z.split()[3].interpolate(Constant(10)) + (u_guess, p_guess) = g.subfunctions + z.subfunctions[0].interpolate(Constant(gamma)) + z.subfunctions[1].assign(u_guess) + z.subfunctions[2].assign(p_guess) + z.subfunctions[3].interpolate(Constant(10)) solver_params = { "snes_monitor": None, diff --git a/tests/regression/test_fieldsplit_split_reorder_bcs.py b/tests/regression/test_fieldsplit_split_reorder_bcs.py index dd6d6f6624..d3ebb06c94 100644 --- a/tests/regression/test_fieldsplit_split_reorder_bcs.py +++ b/tests/regression/test_fieldsplit_split_reorder_bcs.py @@ -118,7 +118,7 @@ def run(solver, solution, permute): solver.solve() sol = solver._problem.u errors = tuple(errornorm(expect, actual, 'L2') for - expect, actual in zip(solution, permute(*sol.split()))) + expect, actual in zip(solution, permute(*sol.subfunctions))) diff = numpy.abs(errors - numpy.asarray([0.02551217479, 0.01991075140, 0.22550499155, diff --git a/tests/regression/test_function_spaces.py b/tests/regression/test_function_spaces.py index 1d33c0f0d6..c4a5e65cc1 100644 --- a/tests/regression/test_function_spaces.py +++ b/tests/regression/test_function_spaces.py @@ -80,14 +80,14 @@ def test_function_space_vector_function_space_differ(mesh): def test_indexed_function_space_index(fs): assert [s.index for s in fs] == list(range(2)) # Create another mixed space in reverse order - fs0, fs1 = fs.split() + fs0, fs1 = fs.subfunctions assert [s.index for s in (fs1 * fs0)] == list(range(2)) # Verify the indices of the original IndexedFunctionSpaces haven't changed assert fs0.index == 0 and fs1.index == 1 def test_mixed_function_space_split(fs): - assert fs.split() == tuple(fs) + assert fs.subfunctions == tuple(fs) def test_function_space_collapse(cg1): diff --git a/tests/regression/test_helmholtz_sphere.py b/tests/regression/test_helmholtz_sphere.py index 0724cd1bbc..403db5afee 100644 --- a/tests/regression/test_helmholtz_sphere.py +++ b/tests/regression/test_helmholtz_sphere.py @@ -54,7 +54,7 @@ def run_helmholtz_mixed_sphere(MeshClass, r, meshd, eltd): 'fieldsplit_0_ksp_type': 'cg', 'fieldsplit_1_ksp_type': 'cg'}) - _, u = soln.split() + _, u = soln.subfunctions f.interpolate(x[0]*x[1]*x[2]/13.0) return errornorm(f, u, degree_rise=0) diff --git a/tests/regression/test_l2pullbacks.py b/tests/regression/test_l2pullbacks.py index a558592e80..4156e46ecc 100644 --- a/tests/regression/test_l2pullbacks.py +++ b/tests/regression/test_l2pullbacks.py @@ -88,7 +88,7 @@ def helmholtz_mixed(r, meshtype, family, hdegree, vdegree=None, meshd=None, usea f.project(sin(x[0]*pi*2)*sin(x[1]*pi*2)) return sqrt(assemble(dot(sol[2] - f, sol[2] - f) * dx)) elif meshtype in ['spherical-quad', 'spherical-tri']: - _, u = sol.split() + _, u = sol.subfunctions f.project(x[0]*x[1]*x[2]/13.0) return errornorm(f, u, degree_rise=0) diff --git a/tests/regression/test_manifolds.py b/tests/regression/test_manifolds.py index a9b69888d3..74fff831ec 100644 --- a/tests/regression/test_manifolds.py +++ b/tests/regression/test_manifolds.py @@ -43,7 +43,7 @@ def run_no_manifold(): solve(a == L, up, bcs=bc, nullspace=nullspace, solver_parameters=params) exact = Function(V1).interpolate(x[0] - 0.5) - u, p = up.split() + u, p = up.subfunctions assert errornorm(exact, p, degree_rise=0) < 1e-8 @@ -85,7 +85,7 @@ def run_manifold(): solve(a == L, up, bcs=bc, nullspace=nullspace, solver_parameters=params) exact = Function(V1).interpolate(x_n[0] - 0.5) - u, p = up.split() + u, p = up.subfunctions assert errornorm(exact, p, degree_rise=0) < 1e-8 diff --git a/tests/regression/test_mixed_interior_facets.py b/tests/regression/test_mixed_interior_facets.py index 4c3515937b..ccc554f9a5 100644 --- a/tests/regression/test_mixed_interior_facets.py +++ b/tests/regression/test_mixed_interior_facets.py @@ -35,7 +35,7 @@ def test_mfs(mesh2D): W = V3 * V1 * V2 u = Function(W) - u0, u1, u2 = u.split() + u0, u1, u2 = u.subfunctions u0.interpolate(Constant(1)) u1.project(Constant((-1.0, -1.0))) u2.interpolate(Constant(1)) diff --git a/tests/regression/test_mixed_tensor.py b/tests/regression/test_mixed_tensor.py index a75e5b2f98..7d4e061e24 100644 --- a/tests/regression/test_mixed_tensor.py +++ b/tests/regression/test_mixed_tensor.py @@ -21,7 +21,7 @@ def test_mass_mixed_tensor(W): a = (inner(u, v) + inner(p, q) + inner(s, t))*dx - V, Q, T = W.split() + V, Q, T = W.subfunctions u = TrialFunction(V) v = TestFunction(V) diff --git a/tests/regression/test_moore_spence.py b/tests/regression/test_moore_spence.py index 482dde83c6..c09223c62f 100644 --- a/tests/regression/test_moore_spence.py +++ b/tests/regression/test_moore_spence.py @@ -58,9 +58,9 @@ def residual(theta, lmbda, ttheta): # Set initial guesses for state, parameter, null eigenmode z = Function(Z) - z.split()[0].assign(th) - z.split()[1].assign(lm) - z.split()[2].assign(eigenmode) + z.subfunctions[0].assign(th) + z.subfunctions[1].assign(lm) + z.subfunctions[2].assign(eigenmode) # Write Moore-Spence system of equations theta, lmbda, phi = split(z) diff --git a/tests/regression/test_mtw.py b/tests/regression/test_mtw.py index 817d583a79..05492ed89d 100755 --- a/tests/regression/test_mtw.py +++ b/tests/regression/test_mtw.py @@ -61,7 +61,7 @@ def test_mtw(): solve(F == 0, up, solver_parameters=params) - u, p = up.split() + u, p = up.subfunctions l2_u.append(errornorm(uex, u)) l2_p.append(errornorm(pex, p)) diff --git a/tests/regression/test_nonlinear_stokes_hdiv.py b/tests/regression/test_nonlinear_stokes_hdiv.py index 5225139e3e..742fa314bb 100644 --- a/tests/regression/test_nonlinear_stokes_hdiv.py +++ b/tests/regression/test_nonlinear_stokes_hdiv.py @@ -92,7 +92,7 @@ def N(u): solve(F == 0, w, bcs=bcs, solver_parameters=params) - u, p = w.split() + u, p = w.subfunctions # test for penetration on bottom assert sqrt(assemble(dot(u, n)**2*ds(3))) < 1e-14 diff --git a/tests/regression/test_nullspace.py b/tests/regression/test_nullspace.py index f962f394b1..35109bd7a7 100644 --- a/tests/regression/test_nullspace.py +++ b/tests/regression/test_nullspace.py @@ -126,7 +126,7 @@ def test_nullspace_mixed(): exact = Function(DG) exact.interpolate(x[1] - 0.5) - sigma, u = w.split() + sigma, u = w.subfunctions assert sqrt(assemble(inner((u - exact), (u - exact))*dx)) < 1e-7 # Now using a Schur complement @@ -141,7 +141,7 @@ def test_nullspace_mixed(): 'fieldsplit_1_ksp_type': 'cg', 'fieldsplit_1_pc_type': 'none'}) - sigma, u = w.split() + sigma, u = w.subfunctions assert sqrt(assemble(inner((u - exact), (u - exact))*dx)) < 5e-8 @@ -313,13 +313,13 @@ def test_near_nullspace_mixed(): bcs = [DirichletBC(W[0].sub(0), 0, (1, 2)), DirichletBC(W[0].sub(1), 0, (3, 4))] rotW = Function(W) - rotV, _ = rotW.split() + rotV, _ = rotW.subfunctions rotV.interpolate(as_vector((-y, x))) c0 = Function(W) - c0V, _ = c0.split() + c0V, _ = c0.subfunctions c1 = Function(W) - c1V, _ = c1.split() + c1V, _ = c1.subfunctions c0V.interpolate(Constant([1., 0.])) c1V.interpolate(Constant([0., 1.])) diff --git a/tests/regression/test_par_loops.py b/tests/regression/test_par_loops.py index 0dba42b00c..9e9a55d5ae 100644 --- a/tests/regression/test_par_loops.py +++ b/tests/regression/test_par_loops.py @@ -66,7 +66,7 @@ def test_mixed_direct_par_loop(f_mixed): """ par_loop((domain, instructions), direct, {'c': (f_mixed, WRITE)}, is_loopy_kernel=True) - assert all(np.allclose(f.dat.data, 1.0) for f in f_mixed.split()) + assert all(np.allclose(f.dat.data, 1.0) for f in f_mixed.subfunctions) @pytest.mark.parametrize('idx', [0, 1]) @@ -123,7 +123,7 @@ def test_indirect_par_loop_read_const_mixed(f_mixed, const): """ par_loop((domain, instructions), dx, {'d': (f_mixed, WRITE), 'constant': (const, READ)}, is_loopy_kernel=True) - assert all(np.allclose(f.dat.data, const.dat.data) for f in f_mixed.split()) + assert all(np.allclose(f.dat.data, const.dat.data) for f in f_mixed.subfunctions) @pytest.mark.parallel(nprocs=2) diff --git a/tests/regression/test_patch_pc.py b/tests/regression/test_patch_pc.py index af9eac1410..cba143ee9e 100644 --- a/tests/regression/test_patch_pc.py +++ b/tests/regression/test_patch_pc.py @@ -53,7 +53,7 @@ def test_jacobi_sor_equivalence(mesh, problem_type, multiplicative): a = (inner(f[i], f[i]) * inner(grad(u), grad(v)))*dx L = inner(Constant(rhs), v)*dx bcs = [DirichletBC(Q, zero(Q.ufl_element().value_shape()), "on_boundary") - for Q in V.split()] + for Q in V.subfunctions] else: a = inner(grad(u), grad(v))*dx L = inner(Constant(rhs), v)*dx diff --git a/tests/regression/test_piola_mixed_fn.py b/tests/regression/test_piola_mixed_fn.py index 58d1c4ffda..743a7ff2e3 100644 --- a/tests/regression/test_piola_mixed_fn.py +++ b/tests/regression/test_piola_mixed_fn.py @@ -34,7 +34,7 @@ def test_sphere_project(): W = U1*U2*U3 f = Function(W) - f1, f2, f3 = f.split() + f1, f2, f3 = f.subfunctions f1.assign(1) f2.assign(2) f3.assign(3) diff --git a/tests/regression/test_point_eval_api.py b/tests/regression/test_point_eval_api.py index 383812ac55..ffe17d6dfc 100644 --- a/tests/regression/test_point_eval_api.py +++ b/tests/regression/test_point_eval_api.py @@ -118,7 +118,7 @@ def test_dont_raise_mixed(): V2 = FunctionSpace(mesh, "RT", 2) V = V1 * V2 f = Function(V) - f1, f2 = f.split() + f1, f2 = f.subfunctions x = SpatialCoordinate(mesh) f1.interpolate(x[0] + 1.2*x[1]) f2.project(as_vector((x[1], 0.8 + x[0]))) diff --git a/tests/regression/test_point_eval_fs.py b/tests/regression/test_point_eval_fs.py index 947a7bc4ce..0ac1a0439a 100644 --- a/tests/regression/test_point_eval_fs.py +++ b/tests/regression/test_point_eval_fs.py @@ -111,7 +111,7 @@ def test_triangle_mixed(mesh_triangle): V2 = FunctionSpace(mesh_triangle, "RT", 2) V = V1 * V2 f = Function(V) - f1, f2 = f.split() + f1, f2 = f.subfunctions x = SpatialCoordinate(mesh_triangle) f1.interpolate(x[0] + 1.2*x[1]) f2.project(as_vector((x[1], 0.8 + x[0]))) diff --git a/tests/regression/test_poisson_mixed_no_bcs.py b/tests/regression/test_poisson_mixed_no_bcs.py index 4246f28515..601fccbfe4 100644 --- a/tests/regression/test_poisson_mixed_no_bcs.py +++ b/tests/regression/test_poisson_mixed_no_bcs.py @@ -55,7 +55,7 @@ def poisson_mixed(size, parameters={}, quadrilateral=False): # Compute solution w = Function(W) solve(a == L, w, solver_parameters=parameters) - sigma, u = w.split() + sigma, u = w.subfunctions # Analytical solution f.interpolate(x[0]*(1-x[0])*x[1]*(1-x[1])) diff --git a/tests/regression/test_poisson_mixed_strong_bcs.py b/tests/regression/test_poisson_mixed_strong_bcs.py index fc46078447..ec4818625d 100644 --- a/tests/regression/test_poisson_mixed_strong_bcs.py +++ b/tests/regression/test_poisson_mixed_strong_bcs.py @@ -55,7 +55,7 @@ def poisson_mixed(size, parameters={}, quadrilateral=False): # Compute solution w = Function(W) solve(a == L, w, bcs=bcs, solver_parameters=parameters) - sigma, u = w.split() + sigma, u = w.subfunctions # Analytical solution f.interpolate(42*x[1]) diff --git a/tests/regression/test_poisson_sphere.py b/tests/regression/test_poisson_sphere.py index 3d768991e6..fba0a7a49c 100644 --- a/tests/regression/test_poisson_sphere.py +++ b/tests/regression/test_poisson_sphere.py @@ -35,7 +35,7 @@ def run_hdiv_l2(MeshClass, refinement, hdiv_space, degree): 'pc_fieldsplit_schur_fact_type': 'FULL', 'fieldsplit_0_ksp_max_it': 100}) - sigma, u = w.split() + sigma, u = w.subfunctions L2_error_u = errornorm(u_exact, u, degree_rise=1) diff --git a/tests/regression/test_projection.py b/tests/regression/test_projection.py index 8efe9f453a..2537b0ae4f 100644 --- a/tests/regression/test_projection.py +++ b/tests/regression/test_projection.py @@ -228,7 +228,7 @@ def test_mixed_projector(mat_type): xs = SpatialCoordinate(m) v = Function(Vc) - v0, v1 = v.split() # sub_functions + v0, v1 = v.subfunctions v0.interpolate(xs[0]*xs[1] + cos(xs[0]+xs[1])) v1.interpolate(xs[0]*xs[1] + sin(xs[0]+xs[1])) mass1 = assemble(sum(split(v))*dx) diff --git a/tests/regression/test_real_space.py b/tests/regression/test_real_space.py index 63e9c50c33..ee5b09a5d3 100644 --- a/tests/regression/test_real_space.py +++ b/tests/regression/test_real_space.py @@ -161,7 +161,7 @@ def poisson(resolution): f = Function(mfs) x = SpatialCoordinate(mesh) - f0, _ = f.split() + f0, _ = f.subfunctions f0.interpolate(cos(x[0])) @@ -194,7 +194,7 @@ def poisson(resolution): f = Function(mfs) x = SpatialCoordinate(mesh) - f0, _ = f.split() + f0, _ = f.subfunctions f0.interpolate(cos(x[0])) @@ -230,7 +230,7 @@ def test_real_space_mixed_assign(): f = Function(W) - q, v = f.split() + q, v = f.subfunctions q.assign(2) g = Function(V) diff --git a/tests/regression/test_split.py b/tests/regression/test_split.py index 4fc9c85754..fcf8f221f8 100644 --- a/tests/regression/test_split.py +++ b/tests/regression/test_split.py @@ -28,7 +28,7 @@ def test_function_split_raises(): f = Function(W) - u, p = f.split() + u, p = f.subfunctions phi = u*dx + p*dx @@ -45,7 +45,7 @@ def test_split_function_derivative(): f = Function(W) - u, p = f.split() + u, p = f.subfunctions f.assign(1) diff --git a/tests/regression/test_stokes_hdiv_parallel.py b/tests/regression/test_stokes_hdiv_parallel.py index 54c3d00d0a..dcace7c6c8 100644 --- a/tests/regression/test_stokes_hdiv_parallel.py +++ b/tests/regression/test_stokes_hdiv_parallel.py @@ -117,7 +117,7 @@ def test_stokes_hdiv_parallel(mat_type, element_pair): solve(a == L, UP, bcs=bcs, nullspace=nullspace, solver_parameters=parameters, appctx=appctx) - u, p = UP.split() + u, p = UP.subfunctions u_error = u - u_exact p_error = p - p_exact err_u.append(sqrt(abs(assemble(inner(u_error, u_error) * dx)))) diff --git a/tests/regression/test_stokes_mini.py b/tests/regression/test_stokes_mini.py index 053787a098..4c0d33bd99 100644 --- a/tests/regression/test_stokes_mini.py +++ b/tests/regression/test_stokes_mini.py @@ -40,7 +40,7 @@ def run_stokes_mini(mat_type, n): w = Function(W) - u, p = w.split() + u, p = w.subfunctions solve(a == L, w, bcs=bcs, solver_parameters={'pc_type': 'fieldsplit', diff --git a/tests/regression/test_vector_laplace_on_quadrilaterals.py b/tests/regression/test_vector_laplace_on_quadrilaterals.py index 602d026a58..bc4875b8ce 100644 --- a/tests/regression/test_vector_laplace_on_quadrilaterals.py +++ b/tests/regression/test_vector_laplace_on_quadrilaterals.py @@ -46,7 +46,7 @@ def vector_laplace(n, degree): 'fieldsplit_1_pc_type': 'lu', 'ksp_monitor': None}) - out_s, out_u = out.split() + out_s, out_u = out.subfunctions return (sqrt(assemble(inner(out_u - exact_u, out_u - exact_u) * dx)), sqrt(assemble(inner(out_s - exact_s, out_s - exact_s) * dx))) diff --git a/tests/regression/test_vfs_component_bcs.py b/tests/regression/test_vfs_component_bcs.py index 13ac2e5684..0748304558 100644 --- a/tests/regression/test_vfs_component_bcs.py +++ b/tests/regression/test_vfs_component_bcs.py @@ -129,7 +129,7 @@ def test_poisson_in_mixed_plus_vfs_components(V, mat_type, make_val): expected = Function(W) - u, p, r = expected.split() + u, p, r = expected.subfunctions x = SpatialCoordinate(V.mesh()) u.interpolate(as_vector((42*x[0], 5*x[1] + 10))) @@ -189,7 +189,7 @@ def test_stokes_component_all(): Ucmp = Function(W) solve(a == L, Ucmp, bcs=bcs_cmp, solver_parameters=params) - for a, b in zip(Uall.split(), Ucmp.split()): + for a, b in zip(Uall.subfunctions, Ucmp.subfunctions): assert np.allclose(a.dat.data_ro, b.dat.data_ro) diff --git a/tests/slate/test_assemble_tensors.py b/tests/slate/test_assemble_tensors.py index 32f2afdced..d972e3588d 100644 --- a/tests/slate/test_assemble_tensors.py +++ b/tests/slate/test_assemble_tensors.py @@ -35,7 +35,7 @@ def function_space(request, mesh): def f(function_space): """Generate a Firedrake function given a particular function space.""" f = Function(function_space) - f_split = f.split() + f_split = f.subfunctions x = SpatialCoordinate(function_space.mesh()) # NOTE: interpolation of UFL expressions into mixed @@ -56,7 +56,7 @@ def f(function_space): def g(function_space): """Generates a Firedrake function given a particular function space.""" g = Function(function_space) - g_split = g.split() + g_split = g.subfunctions x = SpatialCoordinate(function_space.mesh()) # NOTE: interpolation of UFL expressions into mixed @@ -169,7 +169,7 @@ def test_mixed_coefficient_scalar(mesh): V = FunctionSpace(mesh, "DG", 0) W = V * V f = Function(W) - g, h = f.split() + g, h = f.subfunctions f.assign(1) assert np.allclose(assemble(Tensor((g + f[0] + h + f[1])*dx)), 4.0) diff --git a/tests/slate/test_darcy_hybridized_mixed.py b/tests/slate/test_darcy_hybridized_mixed.py index 2e5b659e3c..d87cd23e80 100644 --- a/tests/slate/test_darcy_hybridized_mixed.py +++ b/tests/slate/test_darcy_hybridized_mixed.py @@ -43,7 +43,7 @@ def test_darcy_flow_hybridization(degree, hdiv_family): 'hybridization': {'ksp_type': 'preonly', 'pc_type': 'lu'}} solve(a == L, w, bcs=bcs, solver_parameters=hybrid_params) - sigma_h, u_h = w.split() + sigma_h, u_h = w.subfunctions w2 = Function(W) sc_params = {'mat_type': 'aij', @@ -51,7 +51,7 @@ def test_darcy_flow_hybridization(degree, hdiv_family): 'pc_type': 'lu', 'pc_factor_mat_solver_type': 'mumps'} solve(a == L, w2, bcs=bcs, solver_parameters=sc_params) - nh_sigma, nh_u = w2.split() + nh_sigma, nh_u = w2.subfunctions # Return the L2 error sigma_err = errornorm(sigma_h, nh_sigma) diff --git a/tests/slate/test_hdg_poisson.py b/tests/slate/test_hdg_poisson.py index 9443798739..86a381c358 100644 --- a/tests/slate/test_hdg_poisson.py +++ b/tests/slate/test_hdg_poisson.py @@ -85,7 +85,7 @@ def run_LDG_H_problem(r, degree, quads=False): solver.solve() # Computed flux, scalar, and trace - q_h, u_h, uhat_h = s.split() + q_h, u_h, uhat_h = s.subfunctions scalar_error = errornorm(a_scalar, u_h, norm_type="L2") flux_error = errornorm(a_flux, q_h, norm_type="L2") diff --git a/tests/slate/test_hybrid_poisson_sphere.py b/tests/slate/test_hybrid_poisson_sphere.py index 297dbe0d54..84da99420e 100644 --- a/tests/slate/test_hybrid_poisson_sphere.py +++ b/tests/slate/test_hybrid_poisson_sphere.py @@ -41,7 +41,7 @@ def nullspace_basis(T): appctx = {'trace_nullspace': nullspace_basis} solve(a == L, w, solver_parameters=params, appctx=appctx) - u_h, _ = w.split() + u_h, _ = w.subfunctions error = errornorm(u_exact, u_h) return error diff --git a/tests/slate/test_matrix_free_hybridization.py b/tests/slate/test_matrix_free_hybridization.py index 10a1d99d81..462d3de7a7 100644 --- a/tests/slate/test_matrix_free_hybridization.py +++ b/tests/slate/test_matrix_free_hybridization.py @@ -35,7 +35,7 @@ def test_matrix_free_hybridization(): 'ksp_rtol': 1e-8, 'mat_type': 'matfree'}} solve(a == L, w, bcs=bcs, solver_parameters=matfree_params) - sigma_h, u_h = w.split() + sigma_h, u_h = w.subfunctions w2 = Function(W) aij_params = {'mat_type': 'matfree', @@ -47,7 +47,7 @@ def test_matrix_free_hybridization(): 'ksp_rtol': 1e-8, 'mat_type': 'aij'}} solve(a == L, w2, bcs=bcs, solver_parameters=aij_params) - _sigma, _u = w2.split() + _sigma, _u = w2.subfunctions # Return the L2 error sigma_err = errornorm(sigma_h, _sigma) diff --git a/tests/slate/test_slate_hybridization.py b/tests/slate/test_slate_hybridization.py index 6f14c8145d..157a254db9 100644 --- a/tests/slate/test_slate_hybridization.py +++ b/tests/slate/test_slate_hybridization.py @@ -63,7 +63,7 @@ def setup_poisson_3D(): W = V * U sigma, u = TrialFunctions(W) tau, v = TestFunctions(W) - V, U = W.split() + V, U = W.subfunctions f = Function(U) x, y, z = SpatialCoordinate(mesh) expr = (1+12*pi*pi)*cos(100*pi*x)*cos(100*pi*y)*cos(100*pi*z) @@ -109,7 +109,7 @@ def test_slate_hybridization(degree, hdiv_family, quadrilateral): 'hybridization': {'ksp_type': 'preonly', 'pc_type': 'lu'}} solve(a == L, w, solver_parameters=params) - sigma_h, u_h = w.split() + sigma_h, u_h = w.subfunctions # (Non-hybrid) Need to slam it with preconditioning due to the # system's indefiniteness @@ -122,7 +122,7 @@ def test_slate_hybridization(degree, hdiv_family, quadrilateral): 'pc_fieldsplit_schur_fact_type': 'FULL', 'fieldsplit_0_ksp_type': 'cg', 'fieldsplit_1_ksp_type': 'cg'}) - nh_sigma, nh_u = w2.split() + nh_sigma, nh_u = w2.subfunctions # Return the L2 error sigma_err = errornorm(sigma_h, nh_sigma) @@ -186,7 +186,7 @@ def test_slate_hybridization_nested_schur(): 'preonly_Shat': False, 'jacobi_Shat': False} builder = solver.snes.ksp.pc.getPythonContext().getSchurComplementBuilder() assert options_check(builder, expected), "Some solver options have not ended up in the PC as wanted." - sigma_h, u_h = w.split() + sigma_h, u_h = w.subfunctions w2 = Function(W) solve(a == L, w2, solver_parameters={'ksp_type': 'preonly', @@ -195,7 +195,7 @@ def test_slate_hybridization_nested_schur(): 'pc_python_type': 'firedrake.HybridizationPC', 'hybridization': {'ksp_type': 'preonly', 'pc_type': 'lu'}}) - nh_sigma, nh_u = w2.split() + nh_sigma, nh_u = w2.subfunctions # Return the L2 error sigma_err = errornorm(sigma_h, nh_sigma) @@ -289,7 +289,7 @@ def test_mixed_poisson_approximated_schur(): builder = solver.snes.ksp.pc.getPythonContext().getSchurComplementBuilder() assert options_check(builder, expected), "Some solver options have not ended up in the PC as wanted." - sigma_h, u_h = w.split() + sigma_h, u_h = w.subfunctions # setup second solver w2 = Function(W) @@ -302,7 +302,7 @@ def test_mixed_poisson_approximated_schur(): 'ksp_rtol': 1e-8, 'mat_type': 'matfree'}} solve(a == L, w2, solver_parameters=aij_params) - _sigma, _u = w2.split() + _sigma, _u = w2.subfunctions # Return the L2 error sigma_err = errornorm(sigma_h, _sigma) @@ -354,7 +354,7 @@ def test_slate_hybridization_jacobi_prec_A00(): builder = solver.snes.ksp.pc.getPythonContext().getSchurComplementBuilder() assert options_check(builder, expected), "Some solver options have not ended up in the PC as wanted." - sigma_h, u_h = w.split() + sigma_h, u_h = w.subfunctions # setup second solver w2 = Function(W) @@ -367,7 +367,7 @@ def test_slate_hybridization_jacobi_prec_A00(): 'pc_type': 'none', 'ksp_rtol': 1e-8, 'mat_type': 'matfree'}}) - nh_sigma, nh_u = w2.split() + nh_sigma, nh_u = w2.subfunctions # Return the L2 error sigma_err = errornorm(sigma_h, nh_sigma) @@ -418,7 +418,7 @@ def test_slate_hybridization_jacobi_prec_schur(): builder = solver.snes.ksp.pc.getPythonContext().getSchurComplementBuilder() assert options_check(builder, expected), "Some solver options have not ended up in the PC as wanted." - sigma_h, u_h = w.split() + sigma_h, u_h = w.subfunctions # setup second solver w2 = Function(W) @@ -431,7 +431,7 @@ def test_slate_hybridization_jacobi_prec_schur(): 'pc_type': 'none', 'ksp_rtol': 1e-8, 'mat_type': 'matfree'}}) - nh_sigma, nh_u = w2.split() + nh_sigma, nh_u = w2.subfunctions # Return the L2 error sigma_err = errornorm(sigma_h, nh_sigma) @@ -484,7 +484,7 @@ def test_mixed_poisson_approximated_schur_jacobi_prec(): builder = solver.snes.ksp.pc.getPythonContext().getSchurComplementBuilder() assert options_check(builder, expected), "Some solver options have not ended up in the PC as wanted." - sigma_h, u_h = w.split() + sigma_h, u_h = w.subfunctions # setup second solver w2 = Function(W) @@ -497,7 +497,7 @@ def test_mixed_poisson_approximated_schur_jacobi_prec(): 'ksp_rtol': 1e-8, 'mat_type': 'matfree'}} solve(a == L, w2, solver_parameters=aij_params) - _sigma, _u = w2.split() + _sigma, _u = w2.subfunctions # Return the L2 error sigma_err = errornorm(sigma_h, _sigma) diff --git a/tests/slate/test_slate_hybridization_extr.py b/tests/slate/test_slate_hybridization_extr.py index ab51428ec1..674b159796 100644 --- a/tests/slate/test_slate_hybridization_extr.py +++ b/tests/slate/test_slate_hybridization_extr.py @@ -47,7 +47,7 @@ def test_hybrid_extr_helmholtz(quad): 'pc_type': 'lu', 'pc_factor_mat_solver_type': 'mumps'}} solve(a == L, w, solver_parameters=params) - sigma_h, u_h = w.split() + sigma_h, u_h = w.subfunctions w2 = Function(W) params2 = {'pc_type': 'fieldsplit', @@ -58,7 +58,7 @@ def test_hybrid_extr_helmholtz(quad): 'fieldsplit_0': {'ksp_type': 'cg'}, 'fieldsplit_1': {'ksp_type': 'cg'}} solve(a == L, w2, solver_parameters=params2) - nh_sigma, nh_u = w2.split() + nh_sigma, nh_u = w2.subfunctions sigma_err = errornorm(sigma_h, nh_sigma) u_err = errornorm(u_h, nh_u) diff --git a/tests/slate/test_slate_hybridized_mixed_bcs.py b/tests/slate/test_slate_hybridized_mixed_bcs.py index 1ee7e8aef5..7f36048531 100644 --- a/tests/slate/test_slate_hybridized_mixed_bcs.py +++ b/tests/slate/test_slate_hybridized_mixed_bcs.py @@ -35,7 +35,7 @@ def test_slate_hybridized_on_boundary(degree, hdiv_family, quadrilateral): bcs = [DirichletBC(W[0], zero(), "on_boundary")] solve(a == L, w, solver_parameters=params, bcs=bcs) - sigma_h, u_h = w.split() + sigma_h, u_h = w.subfunctions # (Non-hybrid) w2 = Function(W) @@ -43,7 +43,7 @@ def test_slate_hybridized_on_boundary(degree, hdiv_family, quadrilateral): solver_parameters={'pc_type': 'lu', 'mat_type': 'aij', 'ksp_type': 'preonly'}, bcs=bcs) - nh_sigma, nh_u = w2.split() + nh_sigma, nh_u = w2.subfunctions # Return the L2 error sigma_err = errornorm(sigma_h, nh_sigma) @@ -87,7 +87,7 @@ def test_slate_hybridized_extruded_bcs(degree, hdiv_family): bcs = [DirichletBC(W[0], zero(), "top"), DirichletBC(W[0], zero(), "bottom")] solve(a == L, w, solver_parameters=params, bcs=bcs) - sigma_h, u_h = w.split() + sigma_h, u_h = w.subfunctions # (Non-hybrid) w2 = Function(W) @@ -95,7 +95,7 @@ def test_slate_hybridized_extruded_bcs(degree, hdiv_family): solver_parameters={'pc_type': 'lu', 'mat_type': 'aij', 'ksp_type': 'preonly'}, bcs=bcs) - nh_sigma, nh_u = w2.split() + nh_sigma, nh_u = w2.subfunctions # Return the L2 error sigma_err = errornorm(sigma_h, nh_sigma) diff --git a/tests/slate/test_slate_mixed_direct.py b/tests/slate/test_slate_mixed_direct.py index 65b10b3d52..1eaf5378a3 100644 --- a/tests/slate/test_slate_mixed_direct.py +++ b/tests/slate/test_slate_mixed_direct.py @@ -60,7 +60,7 @@ def test_slate_mixed_vector(Wd): expect = assemble(action(a, f)) - for c, e in zip(C.split(), expect.split()): + for c, e in zip(C.subfunctions, expect.subfunctions): assert numpy.allclose(c.dat.data_ro, e.dat.data_ro) diff --git a/tests/slate/test_variational_prb.py b/tests/slate/test_variational_prb.py index 7cdeb0d601..650ddf16a1 100644 --- a/tests/slate/test_variational_prb.py +++ b/tests/slate/test_variational_prb.py @@ -61,7 +61,7 @@ def test_lvp_equiv_hdg(degree, nested, elimination): ref_solver = LinearVariationalSolver(ref_problem, solver_parameters=params) ref_solver.solve() - _, __, uhat_ref = s.split() + _, __, uhat_ref = s.subfunctions # Now using Slate expressions only _O = Tensor(a) diff --git a/tests/vertexonly/test_interpolation_from_parent.py b/tests/vertexonly/test_interpolation_from_parent.py index 0cb13682c4..a0f01dcf85 100644 --- a/tests/vertexonly/test_interpolation_from_parent.py +++ b/tests/vertexonly/test_interpolation_from_parent.py @@ -252,7 +252,7 @@ def test_mixed_function_interpolation(parentmesh, vertexcoords, tfs): W = W1 * W2 x = SpatialCoordinate(parentmesh) v = Function(V) - v1, v2 = v.split() + v1, v2 = v.subfunctions # Get Function in V1 # use outer product to check Regge works expr1 = outer(x, x) @@ -268,7 +268,7 @@ def test_mixed_function_interpolation(parentmesh, vertexcoords, tfs): # Interpolate Function in V into W w_v = interpolate(v, W) # Split result and check - w_v1, w_v2 = w_v.split() + w_v1, w_v2 = w_v.subfunctions assert np.allclose(w_v1.dat.data_ro, result1) assert np.allclose(w_v2.dat.data_ro, result2) # try and make reusable Interpolator from V to W @@ -276,7 +276,7 @@ def test_mixed_function_interpolation(parentmesh, vertexcoords, tfs): w_v = Function(W) A_w.interpolate(v, output=w_v) # Split result and check - w_v1, w_v2 = w_v.split() + w_v1, w_v2 = w_v.subfunctions assert np.allclose(w_v1.dat.data_ro, result1) assert np.allclose(w_v2.dat.data_ro, result2) # Enough tests - don't both using it again for a different Function in V From d84d5d6708f29ca59a76ba70196c92ed7e899f93 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Wed, 8 Feb 2023 16:42:53 +0000 Subject: [PATCH 14/14] Add FacetSplitPC (#2708) * Add FacetSplitPC --- firedrake/preconditioners/__init__.py | 1 + firedrake/preconditioners/facet_split.py | 265 +++++++++++++++++++++++ tests/regression/test_facet_split.py | 63 ++++++ 3 files changed, 329 insertions(+) create mode 100644 firedrake/preconditioners/facet_split.py create mode 100644 tests/regression/test_facet_split.py diff --git a/firedrake/preconditioners/__init__.py b/firedrake/preconditioners/__init__.py index 3f67ed1718..e63a3c38f1 100644 --- a/firedrake/preconditioners/__init__.py +++ b/firedrake/preconditioners/__init__.py @@ -10,3 +10,4 @@ from firedrake.preconditioners.hypre_ams import * # noqa: F401 from firedrake.preconditioners.hypre_ads import * # noqa: F401 from firedrake.preconditioners.fdm import * # noqa: F401 +from firedrake.preconditioners.facet_split import * # noqa: F401 diff --git a/firedrake/preconditioners/facet_split.py b/firedrake/preconditioners/facet_split.py new file mode 100644 index 0000000000..a1d134f206 --- /dev/null +++ b/firedrake/preconditioners/facet_split.py @@ -0,0 +1,265 @@ +from functools import partial +from mpi4py import MPI +from pyop2 import op2, PermutedMap +from firedrake.petsc import PETSc +from firedrake.preconditioners.base import PCBase +import firedrake.dmhooks as dmhooks +import numpy + +__all__ = ['FacetSplitPC'] + + +class FacetSplitPC(PCBase): + """ A preconditioner that splits a function into interior and facet DOFs. + + Internally this creates a PETSc PC object that can be controlled + by options using the extra options prefix ``facet_``. + + This allows for statically-condensed preconditioners to be applied to + linear systems involving the matrix applied to the full set of DOFs. Code + generated for the matrix-free operator evaluation in the space with full + DOFs will run faster than the one with interior-facet decoposition, since + the full element has a simpler structure. + """ + + needs_python_pmat = False + _prefix = "facet_" + + _permutation_cache = {} + + def get_permutation(self, V, W): + key = (V, W) + if key not in self._permutation_cache: + indices = get_permutation_map(V, W) + if V._comm.allreduce(numpy.all(indices[:-1] <= indices[1:]), MPI.PROD): + self._permutation_cache[key] = None + else: + self._permutation_cache[key] = indices + return self._permutation_cache[key] + + def initialize(self, pc): + from ufl import RestrictedElement, MixedElement, TensorElement, VectorElement + from firedrake import FunctionSpace, TestFunctions, TrialFunctions + from firedrake.assemble import allocate_matrix, TwoFormAssembler + + _, P = pc.getOperators() + appctx = self.get_appctx(pc) + fcp = appctx.get("form_compiler_parameters") + + prefix = pc.getOptionsPrefix() + options_prefix = prefix + self._prefix + options = PETSc.Options(options_prefix) + mat_type = options.getString("mat_type", "submatrix") + + if P.getType() == "python": + ctx = P.getPythonContext() + a = ctx.a + bcs = tuple(ctx.row_bcs) + else: + ctx = dmhooks.get_appctx(pc.getDM()) + a = ctx.Jp or ctx.J + bcs = tuple(ctx._problem.bcs) + + V = a.arguments()[-1].function_space() + assert len(V) == 1, "Interior-facet decomposition of mixed elements is not supported" + + def restrict(ele, restriction_domain): + if isinstance(ele, VectorElement): + return type(ele)(restrict(ele._sub_element, restriction_domain), dim=ele.num_sub_elements()) + elif isinstance(ele, TensorElement): + return type(ele)(restrict(ele._sub_element, restriction_domain), shape=ele._shape, symmetry=ele._symmety) + else: + return RestrictedElement(ele, restriction_domain) + + # W = V[interior] * V[facet] + W = FunctionSpace(V.mesh(), MixedElement([restrict(V.ufl_element(), d) for d in ("interior", "facet")])) + assert W.dim() == V.dim(), "Dimensions of the original and decomposed spaces do not match" + + mixed_operator = a(sum(TestFunctions(W)), sum(TrialFunctions(W)), coefficients={}) + mixed_bcs = tuple(bc.reconstruct(V=W[-1], g=0) for bc in bcs) + + self.perm = None + self.iperm = None + indices = self.get_permutation(V, W) + if indices is not None: + self.perm = PETSc.IS().createGeneral(indices, comm=V._comm) + self.iperm = self.perm.invertPermutation() + + if mat_type != "submatrix": + self.mixed_op = allocate_matrix(mixed_operator, + bcs=mixed_bcs, + form_compiler_parameters=fcp, + mat_type=mat_type, + options_prefix=options_prefix) + self._assemble_mixed_op = TwoFormAssembler(mixed_operator, tensor=self.mixed_op, + form_compiler_parameters=fcp, + bcs=mixed_bcs).assemble + self._assemble_mixed_op() + mixed_opmat = self.mixed_op.petscmat + + def _permute_nullspace(nsp): + if not (nsp.handle and self.iperm): + return nsp + vecs = [vec.duplicate() for vec in nsp.getVecs()] + for vec in vecs: + vec.permute(self.iperm) + return PETSc.NullSpace().create(constant=nsp.hasConstant(), vectors=vecs, comm=nsp.getComm()) + + mixed_opmat.setNullSpace(_permute_nullspace(P.getNullSpace())) + mixed_opmat.setNearNullSpace(_permute_nullspace(P.getNearNullSpace())) + mixed_opmat.setTransposeNullSpace(_permute_nullspace(P.getTransposeNullSpace())) + elif self.perm: + self._permute_op = partial(PETSc.Mat().createSubMatrixVirtual, P, self.iperm, self.iperm) + mixed_opmat = self._permute_op() + else: + mixed_opmat = P + + # Internally, we just set up a PC object that the user can configure + # however from the PETSc command line. Since PC allows the user to specify + # a KSP, we can do iterative by -facet_pc_type ksp. + scpc = PETSc.PC().create(comm=pc.comm) + scpc.incrementTabLevel(1, parent=pc) + + # We set a DM and an appropriate SNESContext on the constructed PC so one + # can do e.g. fieldsplit. + mixed_dm = W.dm + self._dm = mixed_dm + + # Create new appctx + self._ctx_ref = self.new_snes_ctx(pc, + mixed_operator, + mixed_bcs, + mat_type, + fcp, + options_prefix=options_prefix) + + scpc.setDM(mixed_dm) + scpc.setOptionsPrefix(options_prefix) + scpc.setOperators(A=mixed_opmat, P=mixed_opmat) + with dmhooks.add_hooks(mixed_dm, self, appctx=self._ctx_ref, save=False): + scpc.setFromOptions() + self.pc = scpc + + def update(self, pc): + if hasattr(self, "mixed_op"): + self._assemble_mixed_op() + elif hasattr(self, "_permute_op"): + for mat in self.pc.getOperators(): + mat.destroy() + P = self._permute_op() + self.pc.setOperators(A=P, P=P) + self.pc.setUp() + + def apply(self, pc, x, y): + if self.perm: + x.permute(self.iperm) + dm = self._dm + with dmhooks.add_hooks(dm, self, appctx=self._ctx_ref): + self.pc.apply(x, y) + if self.perm: + x.permute(self.perm) + y.permute(self.perm) + + def applyTranspose(self, pc, x, y): + if self.perm: + x.permute(self.iperm) + dm = self._dm + with dmhooks.add_hooks(dm, self, appctx=self._ctx_ref): + self.pc.applyTranspose(x, y) + if self.perm: + x.permute(self.perm) + y.permute(self.perm) + + def view(self, pc, viewer=None): + super(FacetSplitPC, self).view(pc, viewer) + if hasattr(self, "pc"): + viewer.printfASCII("PC using interior-facet decomposition\n") + self.pc.view(viewer) + + def destroy(self, pc): + if hasattr(self, "pc"): + if hasattr(self, "_permute_op"): + for mat in self.pc.getOperators(): + mat.destroy() + self.pc.destroy() + if hasattr(self, "iperm"): + if self.iperm: + self.iperm.destroy() + if hasattr(self, "perm"): + if self.perm: + self.perm.destroy() + + +def split_dofs(elem): + """ Split DOFs into interior and facet DOF, where facets are sorted by entity. + """ + entity_dofs = elem.entity_dofs() + ndim = elem.cell.get_spatial_dimension() + edofs = [[], []] + for key in sorted(entity_dofs.keys()): + vals = entity_dofs[key] + edim = key + try: + edim = sum(edim) + except TypeError: + pass + for k in sorted(vals.keys()): + edofs[edim < ndim].extend(sorted(vals[k])) + return tuple(numpy.array(e, dtype=PETSc.IntType) for e in edofs) + + +def restricted_dofs(celem, felem): + """ Find which DOFs from felem are on celem + :arg celem: the restricted :class:`finat.FiniteElement` + :arg felem: the unrestricted :class:`finat.FiniteElement` + :returns: :class:`numpy.array` with indices of felem that correspond to celem + """ + csplit = split_dofs(celem) + fsplit = split_dofs(felem) + if len(csplit[0]) and len(csplit[1]): + csplit = [numpy.concatenate(csplit)] + fsplit = [numpy.concatenate(fsplit)] + + k = len(csplit[0]) == 0 + if len(csplit[k]) != len(fsplit[k]): + raise ValueError("Finite elements have different DOFs") + perm = numpy.empty_like(csplit[k]) + perm[csplit[k]] = numpy.arange(len(perm), dtype=perm.dtype) + return fsplit[k][perm] + + +def get_permutation_map(V, W): + perm = numpy.empty((V.dof_count, ), dtype=PETSc.IntType) + perm.fill(-1) + vdat = V.make_dat(val=perm) + + offset = 0 + wdats = [] + for Wsub in W: + val = numpy.arange(offset, offset + Wsub.dof_count, dtype=PETSc.IntType) + wdats.append(Wsub.make_dat(val=val)) + offset += Wsub.dof_dset.layout_vec.sizes[0] + + sizes = [Wsub.finat_element.space_dimension() * Wsub.value_size for Wsub in W] + eperm = numpy.concatenate([restricted_dofs(Wsub.finat_element, V.finat_element) for Wsub in W]) + pmap = PermutedMap(V.cell_node_map(), eperm) + + kernel_code = f""" + void permutation(PetscInt *restrict x, + const PetscInt *restrict xi, + const PetscInt *restrict xf){{ + for(PetscInt i=0; i<{sizes[0]}; i++) x[i] = xi[i]; + for(PetscInt i=0; i<{sizes[1]}; i++) x[i+{sizes[0]}] = xf[i]; + return; + }} + """ + kernel = op2.Kernel(kernel_code, "permutation", requires_zeroed_output_arguments=False) + op2.par_loop(kernel, V.mesh().cell_set, + vdat(op2.WRITE, pmap), + wdats[0](op2.READ, W[0].cell_node_map()), + wdats[1](op2.READ, W[1].cell_node_map())) + + own = V.dof_dset.layout_vec.sizes[0] + perm = perm.reshape((-1, )) + perm = V.dof_dset.lgmap.apply(perm, result=perm) + return perm[:own] diff --git a/tests/regression/test_facet_split.py b/tests/regression/test_facet_split.py new file mode 100644 index 0000000000..bcb9e85e2b --- /dev/null +++ b/tests/regression/test_facet_split.py @@ -0,0 +1,63 @@ +import pytest +from firedrake import * + + +@pytest.mark.parametrize(['quadrilateral', 'ptype'], + [(q, p) + for q in [True, False] + for p in ["lu", "jacobi"]]) +def test_facet_split(quadrilateral, ptype): + if ptype == "lu": + parameters = { + "mat_type": "matfree", + "ksp_type": "preonly", + "pc_type": "python", + "pc_python_type": "firedrake.FacetSplitPC", + "facet": { + "mat_type": "aij", + "pc_type": "fieldsplit", + "pc_fieldsplit_type": "schur", + "pc_fieldsplit_schur_fact_type": "full", + "pc_fieldsplit_schur_precondition": "selfp", + "fieldsplit_0_pc_type": "jacobi", + "fieldsplit_1_pc_type": "lu", + }, + } + elif ptype == "jacobi": + parameters = { + "mat_type": "matfree", + "ksp_type": "preonly", + "pc_type": "python", + "pc_python_type": "firedrake.FacetSplitPC", + "facet": { + "mat_type": "submatrix", + "pc_type": "fieldsplit", + "pc_fieldsplit_type": "schur", + "pc_fieldsplit_schur_fact_type": "full", + "pc_fieldsplit_schur_precondition": "a11", + "fieldsplit_0_pc_type": "jacobi", + "fieldsplit_1_pc_type": "jacobi", + "fieldsplit_1_ksp_type": "cg", + "fieldsplit_1_ksp_rtol": 1E-12, + "fieldsplit_1_ksp_atol": 1E-10, + }, + } + + r = 2 + variant = "fdm" if quadrilateral else None + mesh = UnitSquareMesh(2 ** r, 2 ** r, quadrilateral=quadrilateral) + + V = FunctionSpace(mesh, FiniteElement("Lagrange", degree=3, variant=variant)) + u = TrialFunction(V) + v = TestFunction(V) + uh = Function(V) + + a = inner(grad(u), grad(v)) * dx + L = inner(Constant(0), v) * dx + x = SpatialCoordinate(mesh) + u_exact = 42 * x[1] + bcs = [DirichletBC(V, Constant(0), 3), + DirichletBC(V, Constant(42), 4)] + + solve(a == L, uh, bcs=bcs, solver_parameters=parameters) + assert sqrt(assemble(inner(uh - u_exact, uh - u_exact) * dx)) < 1E-10