diff --git a/demos/multicomponent/multicomponent.py.rst b/demos/multicomponent/multicomponent.py.rst index 04178b608c..c21ebbc211 100644 --- a/demos/multicomponent/multicomponent.py.rst +++ b/demos/multicomponent/multicomponent.py.rst @@ -428,7 +428,15 @@ only apply to the normal component of H(div) functions. Elsewhere on the boundary we enforce :math:`J_i \cdot N = 0`. Finally, instead of specifying the value of the barycentric velocity :math:`v` on the inflows and outflows, we enforce :math:`v = \rho^{-1}(G_1 + G_2)`. Boundary conditions that couple -unknowns and/or are nonlinear must be implemented with :class:`~.EquationBC` instead of :class:`~.DirichletBC`. :: +unknowns and/or are nonlinear must be implemented with :class:`~.EquationBC` instead of :class:`~.DirichletBC`. +Since the barycentric velocity is in :math:`[\textrm{CG}_k]^2` it has degrees of freedom located at the points +where the inlets/outlet meet the walls. Even though the boundary conditions on the inlets/outlet +and the walls are compatible at these points, :class:`~.EquationBC` enforces the boundary conditions weakly whilst +:class:`~.DirichletBC` enforces them strongly. Hence, to avoid ambiguity, we set Dirichlet boundary +conditions on the boundary of the boundary, i.e. the points where the inlets/outlet meet the walls. +To enforce these Dirichlet boundary conditions, tuples with the numbers of the boundary edges coincidental +to these points need to be constructed first. This is then passed on to a Dirichlet boundary condition +which is passed on to :class:`~.EquationBC`.:: # Reference species velocities, which we choose to symmetrize so that the molar fluxes agree v_ref_1 = Constant(0.4e-6) # Reference inflow velocity of benzene, m / s @@ -449,7 +457,9 @@ unknowns and/or are nonlinear must be implemented with :class:`~.EquationBC` ins # Boundary conditions on the barycentric velocity are enforced via EquationBC bc_data = {inlet_1_id: rho_v_inflow_1_bc_func, inlet_2_id: rho_v_inflow_2_bc_func, outlet_id: rho_v_outflow_bc_func} F_bc = sum(inner(v - rho_inv * flux, u) * ds(*subdomain) for subdomain, flux in bc_data.items()) - v_bc = EquationBC(F_bc == 0, solution, (*inlet_1_id, *inlet_2_id, *outlet_id), V=Z_h.sub(2)) + ridges = [(i, j) for i in [*inlet_1_id, *inlet_2_id, *outlet_id] for j in walls_ids] + sub_bcs = [DirichletBC(Z_h.sub(2), 0, ridges)] # bc on the boundary of the boundary + v_bc = EquationBC(F_bc == 0, solution, (*inlet_1_id, *inlet_2_id, *outlet_id), V=Z_h.sub(2), bcs=sub_bcs) # The boundary conditions on the fluxes and barycentric velocity # Note that BCs on H(div) spaces only apply to the normal component diff --git a/firedrake/adjoint_utils/function.py b/firedrake/adjoint_utils/function.py index 5d8a211bd7..a4b83b9828 100644 --- a/firedrake/adjoint_utils/function.py +++ b/firedrake/adjoint_utils/function.py @@ -270,7 +270,7 @@ def _ad_dot(self, other, options=None): return assemble(firedrake.inner(self, other)*firedrake.dx) elif riesz_representation == "H1": return assemble((firedrake.inner(self, other) - + firedrake.inner(firedrake.grad(self), other))*firedrake.dx) + + firedrake.inner(firedrake.grad(self), firedrake.grad(other)))*firedrake.dx) else: raise NotImplementedError( "Unknown Riesz representation %s" % riesz_representation) diff --git a/firedrake/assemble.py b/firedrake/assemble.py index d58a4ed6a7..95fa61418c 100644 --- a/firedrake/assemble.py +++ b/firedrake/assemble.py @@ -1255,6 +1255,9 @@ def _apply_bc(self, tensor, bc, u=None): bc.apply(r, u=u) elif isinstance(bc, EquationBCSplit): bc.zero(tensor) + if isinstance(bc.f, ufl.ZeroBaseForm) or bc.f.empty(): + # form is empty, do nothing + return OneFormAssembler(bc.f, bcs=bc.bcs, form_compiler_parameters=self._form_compiler_params, needs_zeroing=False, diff --git a/firedrake/bcs.py b/firedrake/bcs.py index 275744c75e..13f0cd1862 100644 --- a/firedrake/bcs.py +++ b/firedrake/bcs.py @@ -639,9 +639,6 @@ def reconstruct(self, field=None, V=None, subu=None, u=None, row_field=None, col rank = len(self.f.arguments()) splitter = ExtractSubBlock() form = splitter.split(self.f, argument_indices=(row_field, col_field)[:rank]) - if isinstance(form, ufl.ZeroBaseForm) or form.empty(): - # form is empty, do nothing - return if u is not None: form = firedrake.replace(form, {self.u: u}) if action_x is not None: diff --git a/firedrake/fml/form_manipulation_language.py b/firedrake/fml/form_manipulation_language.py index 8eebc94e49..dafba9d99a 100644 --- a/firedrake/fml/form_manipulation_language.py +++ b/firedrake/fml/form_manipulation_language.py @@ -93,7 +93,7 @@ class Term(object): __slots__ = ["form", "labels"] - def __init__(self, form: ufl.Form, label_dict: Mapping = None): + def __init__(self, form: ufl.BaseForm, label_dict: Mapping = None): """ Parameters @@ -216,6 +216,9 @@ def __mul__( __rmul__ = __mul__ + def __neg__(self): + return -1 * self + def __truediv__( self, other: Union[float, Constant, ufl.algebra.Product] @@ -274,7 +277,7 @@ def __init__(self, *terms: Sequence[Term]): def __add__( self, - other: Union[ufl.Form, Term, "LabelledForm"] + other: Union[ufl.BaseForm, Term, "LabelledForm"] ) -> "LabelledForm": """Add a form, term or labelled form to this labelled form. @@ -289,12 +292,12 @@ def __add__( A labelled form containing the terms. """ - if isinstance(other, ufl.Form): - return LabelledForm(*self, Term(other)) - elif type(other) is Term: + if type(other) is Term: return LabelledForm(*self, other) elif type(other) is LabelledForm: return LabelledForm(*self, *other) + elif isinstance(other, ufl.BaseForm): + return LabelledForm(*self, Term(other)) elif other is None: return self else: @@ -304,7 +307,7 @@ def __add__( def __sub__( self, - other: Union[ufl.Form, Term, "LabelledForm"] + other: Union[ufl.BaseForm, Term, "LabelledForm"] ) -> "LabelledForm": """Subtract a form, term or labelled form from this labelled form. @@ -320,15 +323,18 @@ def __sub__( """ if type(other) is Term: - return LabelledForm(*self, Constant(-1.)*other) + return LabelledForm(*self, -other) elif type(other) is LabelledForm: - return LabelledForm(*self, *[Constant(-1.)*t for t in other]) + return LabelledForm(*self, *[-t for t in other]) elif other is None: return self else: # Make new Term for other and subtract it return LabelledForm(*self, Term(Constant(-1.)*other)) + def __rsub__(self, other): + return other + (-self) + def __mul__( self, other: Union[float, Constant, ufl.algebra.Product] @@ -371,6 +377,9 @@ def __truediv__( __rmul__ = __mul__ + def __neg__(self): + return -1 * self + def __iter__(self) -> Sequence: """Iterable of the terms in the labelled form.""" return iter(self.terms) @@ -427,7 +436,7 @@ def label_map( return new_labelled_form @property - def form(self) -> ufl.Form: + def form(self) -> ufl.BaseForm: """Provide the whole form from the labelled form. Raises @@ -437,7 +446,7 @@ def form(self) -> ufl.Form: Returns ------- - ufl.Form + ufl.BaseForm The whole form corresponding to all the terms. """ @@ -479,7 +488,7 @@ def __init__( def __call__( self, - target: Union[ufl.Form, Term, LabelledForm], + target: Union[ufl.BaseForm, Term, LabelledForm], value: Any = None ) -> Union[Term, LabelledForm]: """Apply the label to a form or term. @@ -494,7 +503,7 @@ def __call__( Raises ------ ValueError - If the `target` is not a ufl.Form, Term or + If the `target` is not a ufl.BaseForm, Term or LabelledForm. Returns @@ -514,7 +523,7 @@ def __call__( self.value = self.default_value if isinstance(target, LabelledForm): return LabelledForm(*(self(t, value) for t in target.terms)) - elif isinstance(target, ufl.Form): + elif isinstance(target, ufl.BaseForm): return LabelledForm(Term(target, {self.label: self.value})) elif isinstance(target, Term): new_labels = target.labels.copy() diff --git a/firedrake/interpolation.py b/firedrake/interpolation.py index c8693ce6db..10e7d4600d 100644 --- a/firedrake/interpolation.py +++ b/firedrake/interpolation.py @@ -441,7 +441,11 @@ def __init__(self, expr: Interpolate): else: self.access = op2.WRITE - if self.target_space.ufl_element().mapping() != "identity": + # TODO check V.finat_element.is_lagrange() once https://github.com/firedrakeproject/fiat/pull/200 is released + target_element = self.target_space.ufl_element() + if not ((isinstance(target_element, MixedElement) + and all(sub.mapping() == "identity" for sub in target_element.sub_elements)) + or target_element.mapping() == "identity"): # Identity mapping between reference cell and physical coordinates # implies point evaluation nodes. raise NotImplementedError( @@ -510,7 +514,8 @@ def _get_symbolic_expressions(self) -> tuple[Interpolate, Interpolate]: elif len(shape) == 1: fs_type = partial(VectorFunctionSpace, dim=shape[0]) else: - fs_type = partial(TensorFunctionSpace, shape=shape) + symmetry = self.target_space.ufl_element().symmetry() + fs_type = partial(TensorFunctionSpace, shape=shape, symmetry=symmetry) # Get expression for point evaluation at the dest_node_coords P0DG_vom = fs_type(vom, "DG", 0) diff --git a/firedrake/mg/embedded.py b/firedrake/mg/embedded.py index d38498b2bd..80af65ce64 100644 --- a/firedrake/mg/embedded.py +++ b/firedrake/mg/embedded.py @@ -10,7 +10,7 @@ __all__ = ("TransferManager", ) -native_families = frozenset(["Lagrange", "Discontinuous Lagrange", "Real", "Q", "DQ", "BrokenElement", "Crouzeix-Raviart"]) +native_families = frozenset(["Lagrange", "Discontinuous Lagrange", "Real", "Q", "DQ", "BrokenElement", "Crouzeix-Raviart", "Kong-Mulder-Veldhuizen"]) alfeld_families = frozenset(["Hsieh-Clough-Tocher", "Reduced-Hsieh-Clough-Tocher", "Johnson-Mercier", "Alfeld-Sorokina", "Arnold-Qin", "Reduced-Arnold-Qin", "Christiansen-Hu", "Guzman-Neilan", "Guzman-Neilan Bubble"]) diff --git a/firedrake/scripts/firedrake-zenodo b/firedrake/scripts/firedrake-zenodo index b328237745..f1b536016c 100755 --- a/firedrake/scripts/firedrake-zenodo +++ b/firedrake/scripts/firedrake-zenodo @@ -26,14 +26,12 @@ DESCRIPTIONS = { "ufl": "UFL: the Unified Form Language", "fiat": "FIAT: the Finite Element Automated Tabulator", "petsc": "PETSc: the Portable, Extensible Toolkit for Scientific Computation", - "loopy": "loopy: Transformation-Based Generation of High-Performance CPU/GPU Code", } PYPI_PACKAGE_NAMES = { "firedrake": "firedrake", "ufl": "fenics-ufl", "fiat": "firedrake-fiat", - "loopy": "loopy", "petsc": "petsc4py", } @@ -295,11 +293,11 @@ def zenodo_records(): result = None i = 1 while True: - if i > 20: + if i > 8000//25: raise RuntimeError("More than 8000 uploads on zenodo?") response = requests.get(f"{ZENODO_URL}/records", params={"q": 'owners:19586 OR owners:19587', "all_versions": True, - "size": 400, + "size": 25, "sort": "mostrecent", "page": i}) if response.status_code == 200: diff --git a/firedrake/supermeshing.py b/firedrake/supermeshing.py index 0fe2cbc3f1..af66bbe718 100644 --- a/firedrake/supermeshing.py +++ b/firedrake/supermeshing.py @@ -20,15 +20,18 @@ from pyop2.compilation import load from pyop2.mpi import COMM_SELF from pyop2.utils import get_petsc_dir +from collections import defaultdict __all__ = ["assemble_mixed_mass_matrix", "intersection_finder"] +# TODO replace with KAIJ (we require petsc4py wrappers) class BlockMatrix(object): - def __init__(self, mat, dimension): + def __init__(self, mat, dimension, block_scale=None): self.mat = mat self.dimension = dimension + self.block_scale = block_scale def mult(self, mat, x, y): sizes = self.mat.getSizes() @@ -41,6 +44,8 @@ def mult(self, mat, x, y): xi = PETSc.Vec().createWithArray(xa, size=sizes[1], comm=x.comm) yi = PETSc.Vec().createWithArray(ya, size=sizes[0], comm=y.comm) self.mat.mult(xi, yi) + if self.block_scale is not None: + yi.scale(self.block_scale[i]) y.array[start::stride] = yi.array_r def multTranspose(self, mat, x, y): @@ -54,6 +59,8 @@ def multTranspose(self, mat, x, y): xi = PETSc.Vec().createWithArray(xa, size=sizes[0], comm=x.comm) yi = PETSc.Vec().createWithArray(ya, size=sizes[1], comm=y.comm) self.mat.multTranspose(xi, yi) + if self.block_scale is not None: + yi.scale(self.block_scale[i]) y.array[start::stride] = yi.array_r @@ -68,14 +75,6 @@ def assemble_mixed_mass_matrix(V_A, V_B): if len(V_A) > 1 or len(V_B) > 1: raise NotImplementedError("Sorry, only implemented for non-mixed spaces") - if V_A.ufl_element().mapping() != "identity" or V_B.ufl_element().mapping() != "identity": - msg = """ -Sorry, only implemented for affine maps for now. To do non-affine, we'd need to -import much more of the assembly engine of UFL/TSFC/etc to do the assembly on -each supermesh cell. -""" - raise NotImplementedError(msg) - mesh_A = V_A.mesh() mesh_B = V_B.mesh() @@ -116,15 +115,39 @@ def likely(cell_A): def likely(cell_A): return cell_map[cell_A] - assert V_A.value_size == V_B.value_size - orig_value_size = V_A.value_size - if V_A.value_size > 1: + assert V_A.block_size == V_B.block_size + orig_block_size = V_A.block_size + + # To deal with symmetry, each block of the mass matrix must be rescaled by the multiplicity + if V_A.ufl_element().mapping() == "symmetries": + symmetry = V_A.ufl_element().symmetry() + assert V_B.ufl_element().mapping() == "symmetries" + assert V_B.ufl_element().symmetry() == symmetry + + multiplicity = defaultdict(int) + for idx in numpy.ndindex(V_A.value_shape): + idx = symmetry.get(idx, idx) + multiplicity[idx] += 1 + + block_scale = tuple(scale for idx, scale in multiplicity.items()) + else: + block_scale = None + + if V_A.block_size > 1: V_A = firedrake.FunctionSpace(mesh_A, V_A.ufl_element().sub_elements[0]) - if V_B.value_size > 1: + if V_B.block_size > 1: V_B = firedrake.FunctionSpace(mesh_B, V_B.ufl_element().sub_elements[0]) - assert V_A.value_size == 1 - assert V_B.value_size == 1 + if V_A.ufl_element().mapping() != "identity" or V_B.ufl_element().mapping() != "identity": + msg = """ +Sorry, only implemented for affine maps for now. To do non-affine, we'd need to +import much more of the assembly engine of UFL/TSFC/etc to do the assembly on +each supermesh cell. +""" + raise NotImplementedError(msg) + + assert V_A.block_size == 1 + assert V_B.block_size == 1 preallocator = PETSc.Mat().create(comm=mesh_A.comm) preallocator.setType(PETSc.Mat.Type.PREALLOCATOR) @@ -155,7 +178,7 @@ def likely(cell_A): onnz = numpy.repeat(onnz, cset.cdim) preallocator.destroy() - assert V_A.value_size == V_B.value_size + assert V_A.block_size == V_B.block_size rdim = V_B.dof_dset.cdim cdim = V_A.dof_dset.cdim @@ -445,16 +468,16 @@ def likely(cell_A): lib.restype = ctypes.c_int ammm(V_A, V_B, likely, node_locations_A, node_locations_B, M_SS, ctypes.addressof(lib), mat) - if orig_value_size == 1: + if orig_block_size == 1: return mat else: (lrows, grows), (lcols, gcols) = mat.getSizes() - lrows *= orig_value_size - grows *= orig_value_size - lcols *= orig_value_size - gcols *= orig_value_size + lrows *= orig_block_size + grows *= orig_block_size + lcols *= orig_block_size + gcols *= orig_block_size size = ((lrows, grows), (lcols, gcols)) - context = BlockMatrix(mat, orig_value_size) + context = BlockMatrix(mat, orig_block_size, block_scale=block_scale) blockmat = PETSc.Mat().createPython(size, context=context, comm=mat.comm) blockmat.setUp() return blockmat diff --git a/scripts/firedrake-check b/scripts/firedrake-check index 7aea8faee5..deeb73ca6e 100644 --- a/scripts/firedrake-check +++ b/scripts/firedrake-check @@ -24,6 +24,10 @@ TESTS = { # near nullspace "tests/firedrake/regression/test_nullspace.py::test_near_nullspace", ), + 2: ( + # HDF5/checkpointing + "tests/firedrake/output/test_io_function.py::test_io_function_base[cell_family_degree2]", + ), 3: ( "tests/firedrake/regression/test_dg_advection.py::test_dg_advection_icosahedral_sphere[nprocs=3]", # vertex-only mesh diff --git a/setup.py b/setup.py index 1a9f9361b7..1d5b588355 100644 --- a/setup.py +++ b/setup.py @@ -238,6 +238,7 @@ def extensions(): "tests/firedrake/regression/test_nullspace.py", "tests/firedrake/regression/test_dg_advection.py", "tests/firedrake/regression/test_interpolate_cross_mesh.py", + "tests/firedrake/output/test_io_function.py", ) diff --git a/tests/firedrake/adjoint/test_reduced_functional.py b/tests/firedrake/adjoint/test_reduced_functional.py index 4e33915258..b749f8fb84 100644 --- a/tests/firedrake/adjoint/test_reduced_functional.py +++ b/tests/firedrake/adjoint/test_reduced_functional.py @@ -282,3 +282,31 @@ def test_real_space_parallel(): Jhat = ReducedFunctional(J, Control(m)) opt = minimize(Jhat) parallel_assert(np.allclose(opt.dat.data_ro, 1)) + + +@pytest.mark.parametrize("riesz_representation", ["l2", "L2", "H1"]) +@pytest.mark.skipcomplex +def test_ad_dot(riesz_representation): + mesh = IntervalMesh(10, 0, 1) + V = FunctionSpace(mesh, "Lagrange", 1) + + c = Constant(1) + f = Function(V) + x = SpatialCoordinate(mesh) + f.interpolate(x[0]) + + u = Function(V) + v = TestFunction(V) + bc = DirichletBC(V, Constant(1), "on_boundary") + + F = inner(grad(u), grad(v))*dx - f**2*v*dx + solve(F == 0, u, bc) + + J = assemble(c**2*u*dx) + Jhat = ReducedFunctional(J, Control(f, riesz_map=riesz_representation)) + dJhat = Jhat.derivative(apply_riesz=True) + + h = Function(V) + h.dat.data[:] = np.random.rand(V.dof_dset.size) + dJdh = dJhat._ad_dot(h, options={'riesz_representation': riesz_representation}) + assert taylor_test(Jhat, f, h, dJdm=dJdh) > 1.9 diff --git a/tests/firedrake/equation_bcs/test_equation_bcs.py b/tests/firedrake/equation_bcs/test_equation_bcs.py index d487aa1e27..53daf97a54 100644 --- a/tests/firedrake/equation_bcs/test_equation_bcs.py +++ b/tests/firedrake/equation_bcs/test_equation_bcs.py @@ -174,8 +174,8 @@ def linear_poisson_mixed(solver_parameters, mesh_num, porder): u1 = cos(2 * pi * y) / 2 n = FacetNormal(mesh) - a = (inner(sigma, tau) + inner(u, div(tau)) + inner(div(sigma), v)) * dx - L = inner(u1, dot(tau, n)) * ds(1) + inner(f, v) * dx + a = inner(sigma, tau)*dx + inner(u, div(tau))*dx + inner(div(sigma), v + div(tau))*dx + L = inner(u1, dot(tau, n)) * ds(1) + inner(f, v + div(tau)) * dx g = 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)]) @@ -331,24 +331,24 @@ def test_EquationBC_mixedpoisson_matfree_fieldsplit(): eq_type = "linear" porder = 0 # Mixed poisson with EquationBCs - # matfree with fieldsplit pc + # matfree with fieldsplit schur direct solver solver_parameters = {'mat_type': mat_type, - 'ksp_type': 'fgmres', - 'ksp_atol': 1e-11, - 'ksp_max_it': 200, + 'ksp_type': 'preonly', 'pc_type': 'fieldsplit', 'pc_fieldsplit_type': 'schur', 'pc_fieldsplit_schur_fact_type': 'full', - 'fieldsplit_0_ksp_type': 'cg', + 'fieldsplit_0_ksp_type': 'preonly', 'fieldsplit_0_pc_type': 'python', 'fieldsplit_0_pc_python_type': 'firedrake.AssembledPC', 'fieldsplit_0_assembled_pc_type': 'lu', - 'fieldsplit_1_ksp_type': 'cg', - 'fieldsplit_1_pc_use_amat': True, + 'fieldsplit_0_assembled_pc_factor_mat_solver_type': 'mumps', + 'fieldsplit_1_ksp_type': 'gmres', 'fieldsplit_1_pc_type': 'python', 'fieldsplit_1_pc_python_type': 'firedrake.MassInvPC', - 'fieldsplit_1_Mp_pc_type': 'icc'} + 'fieldsplit_1_Mp_pc_type': 'lu', + 'fieldsplit_1_Mp_pc_factor_mat_solver_type': 'mumps', + } err = [] mesh_sizes = [16, 32] if eq_type == "linear": diff --git a/tests/firedrake/multigrid/test_grid_transfer.py b/tests/firedrake/multigrid/test_grid_transfer.py index f4c689a1c6..86af28b455 100644 --- a/tests/firedrake/multigrid/test_grid_transfer.py +++ b/tests/firedrake/multigrid/test_grid_transfer.py @@ -66,10 +66,7 @@ def shape(request): @pytest.fixture(params=["injection", "restriction", "prolongation"]) def transfer_type(request, hierarchy): - if not hierarchy.nested and request.param == "injection": - return pytest.mark.xfail(reason="Supermesh projections not implemented yet")(request.param) - else: - return request.param + return request.param @pytest.fixture @@ -167,19 +164,35 @@ def functional(victim, dual): assert numpy.allclose(fine_functional, coarse_functional) -def test_grid_transfer(hierarchy, shape, space, degrees, transfer_type): - if not hierarchy.nested and transfer_type == "injection": +def run_transfer(mh, shp, family, deg, transfer_op): + if not mh.nested and mh.refinements_per_level > 1: pytest.skip("Not implemented") - if transfer_type == "injection": - if space in {"DG", "DQ"} and complex_mode: + if transfer_op == "injection": + if not mh.nested: + pytest.skip("Supermesh projections not implemented yet") + if family in {"DG", "DQ"} and complex_mode: with pytest.raises(NotImplementedError): - run_injection(hierarchy, shape, space, degrees) + run_injection(mh, shp, family, deg) else: - run_injection(hierarchy, shape, space, degrees) - elif transfer_type == "restriction": - run_restriction(hierarchy, shape, space, degrees) - elif transfer_type == "prolongation": - run_prolongation(hierarchy, shape, space, degrees) + run_injection(mh, shp, family, deg) + elif transfer_op == "restriction": + run_restriction(mh, shp, family, deg) + elif transfer_op == "prolongation": + run_prolongation(mh, shp, family, deg) + else: + raise ValueError(f"Invalid transfer {transfer_op}") + + +def test_grid_transfer(hierarchy, shape, space, degrees, transfer_type): + run_transfer(hierarchy, shape, space, degrees, transfer_type) + + +@pytest.mark.parallel(nprocs=2) +def test_grid_transfer_parallel(hierarchy, transfer_type): + space = "CG" + degrees = (1, 2, 3) + shape = "scalar" + run_transfer(hierarchy, shape, space, degrees, transfer_type) @pytest.mark.parallel([1, 2]) @@ -191,35 +204,18 @@ def test_grid_transfer_symmetric(transfer_type): space = "Lagrange" degrees = (1,) shape = "symmetric-tensor" - if transfer_type == "injection": - if space in {"DG", "DQ"} and complex_mode: - with pytest.raises(NotImplementedError): - run_injection(hierarchy, shape, space, degrees) - else: - run_injection(hierarchy, shape, space, degrees) - elif transfer_type == "restriction": - run_restriction(hierarchy, shape, space, degrees) - elif transfer_type == "prolongation": - run_prolongation(hierarchy, shape, space, degrees) + run_transfer(hierarchy, shape, space, degrees, transfer_type) -@pytest.mark.parallel(nprocs=2) -def test_grid_transfer_parallel(hierarchy, transfer_type): - space = "CG" - degrees = (1, 2, 3) +@pytest.mark.parametrize("transfer_type", ["prolongation", "restriction", "injection"]) +def test_grid_transfer_KMV(transfer_type): + base = UnitSquareMesh(3, 3) + hierarchy = MeshHierarchy(base, 1) + + space = "KMV" + degrees = (2, 5) shape = "scalar" - if not hierarchy.nested and hierarchy.refinements_per_level > 1: - pytest.skip("Not implemented") - if transfer_type == "injection": - if space in {"DG", "DQ"} and complex_mode: - with pytest.raises(NotImplementedError): - run_injection(hierarchy, shape, space, degrees) - else: - run_injection(hierarchy, shape, space, degrees) - elif transfer_type == "restriction": - run_restriction(hierarchy, shape, space, degrees) - elif transfer_type == "prolongation": - run_prolongation(hierarchy, shape, space, degrees) + run_transfer(hierarchy, shape, space, degrees, transfer_type) @pytest.fixture(params=["interval-interval", @@ -283,16 +279,7 @@ def test_grid_transfer_deformed(deformed_hierarchy, deformed_transfer_type): shape = "scalar" if not deformed_hierarchy.nested and deformed_transfer_type == "injection": pytest.skip("Not implemented") - if deformed_transfer_type == "injection": - if space in {"DG", "DQ"} and complex_mode: - with pytest.raises(NotImplementedError): - run_injection(deformed_hierarchy, shape, space, degrees[:1]) - else: - run_injection(deformed_hierarchy, shape, space, degrees[:1]) - elif deformed_transfer_type == "restriction": - run_restriction(deformed_hierarchy, shape, space, degrees) - elif deformed_transfer_type == "prolongation": - run_prolongation(deformed_hierarchy, shape, space, degrees) + run_transfer(deformed_hierarchy, shape, space, degrees[:1], deformed_transfer_type) @pytest.fixture(params=["interval", "triangle", "quadrilateral", "tetrahedron"], scope="module") diff --git a/tests/firedrake/output/test_hdf5file_checkpoint.py b/tests/firedrake/output/test_hdf5file_checkpoint.py index 22bfc453cd..5d8cf7436d 100644 --- a/tests/firedrake/output/test_hdf5file_checkpoint.py +++ b/tests/firedrake/output/test_hdf5file_checkpoint.py @@ -38,7 +38,8 @@ def f(): return f -def run_write_read(mesh, fs, degree, dumpfile): +@pytest.mark.parallel([1, 2]) +def test_write_read(mesh, fs, degree, dumpfile): V = FunctionSpace(mesh, fs, degree) @@ -79,15 +80,6 @@ def test_checkpoint_fails_for_non_function(dumpfile): h5.write(np.arange(10), "/solution") -def test_write_read(mesh, fs, degree, dumpfile): - run_write_read(mesh, fs, degree, dumpfile) - - -@pytest.mark.parallel(nprocs=2) -def test_write_read_parallel(mesh, fs, degree, dumpfile): - run_write_read(mesh, fs, degree, dumpfile) - - def test_checkpoint_read_not_exist_ioerror(dumpfile): with pytest.raises(IOError): with HDF5File(dumpfile, file_mode="r"): diff --git a/tests/firedrake/output/test_io_function.py b/tests/firedrake/output/test_io_function.py index 7651011759..0ea6812428 100644 --- a/tests/firedrake/output/test_io_function.py +++ b/tests/firedrake/output/test_io_function.py @@ -148,7 +148,7 @@ def _load_check_save_functions(filename, func_name, comm, method, mesh_name, var afile.save_function(fB) -@pytest.mark.parallel(nprocs=2) +@pytest.mark.parallel(2) @pytest.mark.parametrize('cell_family_degree', [ ("triangle_small", "P", 1), ("triangle_small", "P", 6), diff --git a/tests/firedrake/regression/test_fml.py b/tests/firedrake/regression/test_fml.py index f639c8e870..f6c02d55a1 100644 --- a/tests/firedrake/regression/test_fml.py +++ b/tests/firedrake/regression/test_fml.py @@ -9,9 +9,12 @@ PeriodicUnitSquareMesh, FunctionSpace, Constant, MixedFunctionSpace, TestFunctions, Function, split, inner, dx, SpatialCoordinate, as_vector, pi, sin, div, - NonlinearVariationalProblem, NonlinearVariationalSolver + NonlinearVariationalProblem, NonlinearVariationalSolver, + UnitIntervalMesh, Cofunction, TestFunction, interpolate ) -from firedrake.fml import subject, replace_subject, keep, drop, Label +from firedrake.fml import subject, replace_subject, keep, drop, Label, LabelledForm + +import pytest def test_fml(): @@ -121,3 +124,44 @@ def test_fml(): X.assign(X_np1) implicit_solver.solve() X.assign(X_np1) + + +@pytest.fixture +def V(): + mesh = UnitIntervalMesh(2) + return FunctionSpace(mesh, "CG", 1) + + +@pytest.mark.parametrize("baseform", ("form", "cofunction", "interpolate")) +def test_fml_baseform(V, baseform): + if baseform == "form": + form = inner(1, TestFunction(V))*dx + + elif baseform == "cofunction": + form = Cofunction(V.dual()) + form.assign(42) + + elif baseform == "interpolate": + W = V.reconstruct(degree=2) + w = Cofunction(W.dual()) + w.assign(42) + form = interpolate(TestFunction(V), w) + + # Test that we can wrap this BaseForm with a Label + has_labels = lambda term: len(term.labels) > 0 + label = Label("label") + F = label(form) + assert isinstance(F, LabelledForm) + assert F.label_map(has_labels, map_if_true=keep, map_if_false=drop).form == form + + # Test that we can linearly combine this BaseForm with a LabelledForm + Fcombs = [form + label(form), + label(form) + form, + label(form) - form] + for F2 in Fcombs: + assert isinstance(F2, LabelledForm) + assert F2.label_map(has_labels, map_if_true=keep, map_if_false=drop).form == form + + F3 = form - label(form) + assert isinstance(F3, LabelledForm) + assert F3.label_map(has_labels, map_if_true=keep, map_if_false=drop).form == -form diff --git a/tests/firedrake/regression/test_interpolate_cross_mesh.py b/tests/firedrake/regression/test_interpolate_cross_mesh.py index a669618ed3..32179172f7 100644 --- a/tests/firedrake/regression/test_interpolate_cross_mesh.py +++ b/tests/firedrake/regression/test_interpolate_cross_mesh.py @@ -406,11 +406,12 @@ def test_exact_refinement(): ) -def test_interpolate_unitsquare_tfs_shape(): +@pytest.mark.parametrize("shape,symmetry", [((1, 2, 3), None), ((3, 3), True)]) +def test_interpolate_unitsquare_tfs_shape(shape, symmetry): m_src = UnitSquareMesh(2, 3) m_dest = UnitSquareMesh(3, 5, quadrilateral=True) - V_src = TensorFunctionSpace(m_src, "CG", 3, shape=(1, 2, 3)) - V_dest = TensorFunctionSpace(m_dest, "CG", 4, shape=(1, 2, 3)) + V_src = TensorFunctionSpace(m_src, "CG", 3, shape=shape, symmetry=symmetry) + V_dest = TensorFunctionSpace(m_dest, "CG", 4, shape=shape, symmetry=symmetry) f_src = Function(V_src) assemble(interpolate(f_src, V_dest)) diff --git a/tests/firedrake/supermesh/test_galerkin_projection.py b/tests/firedrake/supermesh/test_galerkin_projection.py index 2770bfa7d2..cc20b8b3b9 100644 --- a/tests/firedrake/supermesh/test_galerkin_projection.py +++ b/tests/firedrake/supermesh/test_galerkin_projection.py @@ -2,6 +2,7 @@ from firedrake.petsc import DEFAULT_DIRECT_SOLVER_PARAMETERS from firedrake.supermeshing import * from itertools import product +from functools import partial import numpy import pytest @@ -14,7 +15,7 @@ def mesh(request): return UnitCubeMesh(3, 2, 1) -@pytest.fixture(params=["scalar", "vector", pytest.param("tensor", marks=pytest.mark.skip(reason="Prolongation fails for tensors"))]) +@pytest.fixture(params=["scalar", "vector", "tensor", "symmetric"]) def shapify(request): if request.param == "scalar": return lambda x: x @@ -22,6 +23,8 @@ def shapify(request): return VectorElement elif request.param == "tensor": return TensorElement + elif request.param == "symmetric": + return partial(TensorElement, symmetry=True) else: raise RuntimeError