From e5f0bbf10749ddbb75ad8611664b3f268cc03560 Mon Sep 17 00:00:00 2001 From: Florian Rathgeber Date: Thu, 30 Jan 2014 16:47:20 +0000 Subject: [PATCH 01/27] Migrate FFC interface from PyOP2 --- firedrake/ffc_interface.py | 107 +++++++++++++ firedrake/firedrake_geometry.h | 274 +++++++++++++++++++++++++++++++++ firedrake/solving.py | 5 +- tests/test_ffc_interface.py | 106 +++++++++++++ 4 files changed, 490 insertions(+), 2 deletions(-) create mode 100644 firedrake/ffc_interface.py create mode 100644 firedrake/firedrake_geometry.h create mode 100644 tests/test_ffc_interface.py diff --git a/firedrake/ffc_interface.py b/firedrake/ffc_interface.py new file mode 100644 index 0000000000..6fbfabd61c --- /dev/null +++ b/firedrake/ffc_interface.py @@ -0,0 +1,107 @@ +"""Provides the interface to FFC for compiling a form, and transforms the FFC- +generated code in order to make it suitable for passing to the backends.""" + +from hashlib import md5 +import os +import tempfile + +from ufl import Form +from ufl.algorithms import as_form +from ffc import default_parameters, compile_form as ffc_compile_form +from ffc import constants + +from pyop2.caching import DiskCached +from pyop2.op2 import Kernel +from pyop2.mpi import MPI +from pyop2.ir.ast_base import PreprocessNode, Root + +_form_cache = {} + +ffc_parameters = default_parameters() +ffc_parameters['write_file'] = False +ffc_parameters['format'] = 'pyop2' +ffc_parameters['pyop2-ir'] = True + +# Include an md5 hash of firedrake_geometry.h in the cache key +with open(os.path.join(os.path.dirname(__file__), 'firedrake_geometry.h')) as f: + _firedrake_geometry_md5 = md5(f.read()).hexdigest() + + +def _check_version(): + from pyop2.version import __compatible_ffc_version_info__ as compatible_version, \ + __compatible_ffc_version__ as version + try: + if constants.PYOP2_VERSION_INFO[:2] == compatible_version[:2]: + return + except AttributeError: + pass + raise RuntimeError("Incompatible PyOP2 version %s and FFC PyOP2 version %s." + % (version, getattr(constants, 'PYOP2_VERSION', 'unknown'))) + + +class FFCKernel(DiskCached): + + _cache = {} + _cachedir = os.path.join(tempfile.gettempdir(), + 'firedrake-ffc-kernel-cache-uid%d' % os.getuid()) + + @classmethod + def _cache_key(cls, form, name): + form_data = form.compute_form_data() + return md5(form_data.signature + name + Kernel._backend.__name__ + + _firedrake_geometry_md5 + constants.FFC_VERSION + + constants.PYOP2_VERSION).hexdigest() + + def __init__(self, form, name): + if self._initialized: + return + + incl = PreprocessNode('#include "firedrake_geometry.h"\n') + inc = [os.path.dirname(__file__)] + forms = ffc_compile_form(form, prefix=name, parameters=ffc_parameters) + fdict = dict((f.name, f) for f in forms) + + kernels = [] + for ida in form.form_data().preprocessed_form.integrals(): + fname = '%s_%s_integral_0_%s' % (name, ida.domain_type(), ida.domain_id()) + # Set optimization options + opts = {} if ida.domain_type() not in ['cell'] else \ + {'licm': False, + 'tile': None, + 'vect': None, + 'ap': False, + 'split': None} + kernels.append(Kernel(Root([incl, fdict[fname]]), fname, opts, inc)) + self.kernels = tuple(kernels) + + self._initialized = True + + +def compile_form(form, name): + """Compile a form using FFC and return a tuple of + :class:`Kernels `.""" + + # Check that we get a Form + if not isinstance(form, Form): + form = as_form(form) + + return FFCKernel(form, name).kernels + + +def clear_cache(): + """Clear the PyOP2 FFC kernel cache.""" + if MPI.comm.rank != 0: + return + if os.path.exists(FFCKernel._cachedir): + import shutil + shutil.rmtree(FFCKernel._cachedir, ignore_errors=True) + _ensure_cachedir() + + +def _ensure_cachedir(): + """Ensure that the FFC kernel cache directory exists.""" + if not os.path.exists(FFCKernel._cachedir) and MPI.comm.rank == 0: + os.makedirs(FFCKernel._cachedir) + +_check_version() +_ensure_cachedir() diff --git a/firedrake/firedrake_geometry.h b/firedrake/firedrake_geometry.h new file mode 100644 index 0000000000..5ef324927f --- /dev/null +++ b/firedrake/firedrake_geometry.h @@ -0,0 +1,274 @@ +/* --- Computation of Jacobian matrices --- */ + +/* compute Jacobian J for interval embedded in R^1 */ +#define compute_jacobian_interval_1d(J, vertex_coordinates) \ + J[0] = vertex_coordinates[1][0] - vertex_coordinates[0][0]; + +/* Compute Jacobian J for interval embedded in R^2 */ +#define compute_jacobian_interval_2d(J, vertex_coordinates) \ + J[0] = vertex_coordinates[1][0] - vertex_coordinates[0][0]; \ + J[1] = vertex_coordinates[3][0] - vertex_coordinates[2][0]; + +/* Compute Jacobian J for quad embedded in R^2 */ +#define compute_jacobian_quad_2d(J, vertex_coordinates) \ + J[0] = vertex_coordinates[2][0] - vertex_coordinates[0][0]; \ + J[1] = vertex_coordinates[1][0] - vertex_coordinates[0][0]; \ + J[2] = vertex_coordinates[6][0] - vertex_coordinates[4][0]; \ + J[3] = vertex_coordinates[5][0] - vertex_coordinates[4][0]; + +/* Compute Jacobian J for quad embedded in R^3 */ +#define compute_jacobian_quad_3d(J, vertex_coordinates) \ + J[0] = vertex_coordinates[2] [0] - vertex_coordinates[0][0]; \ + J[1] = vertex_coordinates[1] [0] - vertex_coordinates[0][0]; \ + J[2] = vertex_coordinates[6] [0] - vertex_coordinates[4][0]; \ + J[3] = vertex_coordinates[5] [0] - vertex_coordinates[4][0]; \ + J[4] = vertex_coordinates[10] [0] - vertex_coordinates[8][0]; \ + J[5] = vertex_coordinates[9][0] - vertex_coordinates[8][0]; + +/* Compute Jacobian J for interval embedded in R^3 */ +#define compute_jacobian_interval_3d(J, vertex_coordinates) \ + J[0] = vertex_coordinates[1][0] - vertex_coordinates[0][0]; \ + J[1] = vertex_coordinates[3][0] - vertex_coordinates[2][0]; \ + J[2] = vertex_coordinates[5][0] - vertex_coordinates[4][0]; + +/* Compute Jacobian J for triangle embedded in R^2 */ +#define compute_jacobian_triangle_2d(J, vertex_coordinates) \ + J[0] = vertex_coordinates[1][0] - vertex_coordinates[0][0]; \ + J[1] = vertex_coordinates[2][0] - vertex_coordinates[0][0]; \ + J[2] = vertex_coordinates[4][0] - vertex_coordinates[3][0]; \ + J[3] = vertex_coordinates[5][0] - vertex_coordinates[3][0]; + +/* Compute Jacobian J for triangle embedded in R^3 */ +#define compute_jacobian_triangle_3d(J, vertex_coordinates) \ + J[0] = vertex_coordinates[1][0] - vertex_coordinates[0][0]; \ + J[1] = vertex_coordinates[2][0] - vertex_coordinates[0][0]; \ + J[2] = vertex_coordinates[4][0] - vertex_coordinates[3][0]; \ + J[3] = vertex_coordinates[5][0] - vertex_coordinates[3][0]; \ + J[4] = vertex_coordinates[7][0] - vertex_coordinates[6][0]; \ + J[5] = vertex_coordinates[8][0] - vertex_coordinates[6][0]; + +/* Compute Jacobian J for tetrahedron embedded in R^3 */ +#define compute_jacobian_tetrahedron_3d(J, vertex_coordinates) \ + J[0] = vertex_coordinates[1] [0] - vertex_coordinates[0][0]; \ + J[1] = vertex_coordinates[2] [0] - vertex_coordinates[0][0]; \ + J[2] = vertex_coordinates[3] [0] - vertex_coordinates[0][0]; \ + J[3] = vertex_coordinates[5] [0] - vertex_coordinates[4][0]; \ + J[4] = vertex_coordinates[6] [0] - vertex_coordinates[4][0]; \ + J[5] = vertex_coordinates[7] [0] - vertex_coordinates[4][0]; \ + J[6] = vertex_coordinates[9] [0] - vertex_coordinates[8][0]; \ + J[7] = vertex_coordinates[10][0] - vertex_coordinates[8][0]; \ + J[8] = vertex_coordinates[11][0] - vertex_coordinates[8][0]; + +/* Compute Jacobian J for tensor product prism embedded in R^3 */ +#define compute_jacobian_prism_3d(J, vertex_coordinates) \ + J[0] = vertex_coordinates[2][0] - vertex_coordinates[0][0]; \ + J[1] = vertex_coordinates[4][0] - vertex_coordinates[0][0]; \ + J[2] = vertex_coordinates[1][0] - vertex_coordinates[0][0]; \ + J[3] = vertex_coordinates[8][0] - vertex_coordinates[6][0]; \ + J[4] = vertex_coordinates[10][0] - vertex_coordinates[6][0]; \ + J[5] = vertex_coordinates[7][0] - vertex_coordinates[6][0]; \ + J[6] = vertex_coordinates[14][0] - vertex_coordinates[12][0]; \ + J[7] = vertex_coordinates[16][0] - vertex_coordinates[12][0]; \ + J[8] = vertex_coordinates[13][0] - vertex_coordinates[12][0]; + +/* Jacobians for interior facets of different sorts */ + +/* Compute Jacobian J for interval embedded in R^1 */ +#define compute_jacobian_interval_int_1d compute_jacobian_interval_1d + +/* Compute Jacobian J for interval embedded in R^2 */ +#define compute_jacobian_interval_int_2d(J, vertex_coordinates) \ + J[0] = vertex_coordinates[1][0] - vertex_coordinates[0][0]; \ + J[1] = vertex_coordinates[5][0] - vertex_coordinates[4][0]; + +/* Compute Jacobian J for quad embedded in R^2 */ +#define compute_jacobian_quad_int_2d(J, vertex_coordinates) \ + J[0] = vertex_coordinates[2] [0] - vertex_coordinates[0][0]; \ + J[1] = vertex_coordinates[1] [0] - vertex_coordinates[0][0]; \ + J[2] = vertex_coordinates[10] [0] - vertex_coordinates[8][0]; \ + J[3] = vertex_coordinates[9][0] - vertex_coordinates[8][0]; + +/* Compute Jacobian J for quad embedded in R^3 */ +#define compute_jacobian_quad_int_3d(J, vertex_coordinates) \ + J[0] = vertex_coordinates[2] [0] - vertex_coordinates[0] [0]; \ + J[1] = vertex_coordinates[1] [0] - vertex_coordinates[0] [0]; \ + J[2] = vertex_coordinates[10] [0] - vertex_coordinates[8] [0]; \ + J[3] = vertex_coordinates[9][0] - vertex_coordinates[8] [0]; \ + J[4] = vertex_coordinates[18][0] - vertex_coordinates[16][0]; \ + J[5] = vertex_coordinates[17][0] - vertex_coordinates[16][0]; + +/* Compute Jacobian J for interval embedded in R^3 */ +#define compute_jacobian_interval_int_3d(J, vertex_coordinates) \ + J[0] = vertex_coordinates[1][0] - vertex_coordinates[0][0]; \ + J[1] = vertex_coordinates[5][0] - vertex_coordinates[4][0]; \ + J[2] = vertex_coordinates[9][0] - vertex_coordinates[8][0]; + +/* Compute Jacobian J for triangle embedded in R^2 */ +#define compute_jacobian_triangle_int_2d(J, vertex_coordinates) \ + J[0] = vertex_coordinates[1][0] - vertex_coordinates[0][0]; \ + J[1] = vertex_coordinates[2][0] - vertex_coordinates[0][0]; \ + J[2] = vertex_coordinates[7][0] - vertex_coordinates[6][0]; \ + J[3] = vertex_coordinates[8][0] - vertex_coordinates[6][0]; + +/* Compute Jacobian J for triangle embedded in R^3 */ +#define compute_jacobian_triangle_int_3d(J, vertex_coordinates) \ + J[0] = vertex_coordinates[1] [0] - vertex_coordinates[0] [0]; \ + J[1] = vertex_coordinates[2] [0] - vertex_coordinates[0] [0]; \ + J[2] = vertex_coordinates[7] [0] - vertex_coordinates[6] [0]; \ + J[3] = vertex_coordinates[8] [0] - vertex_coordinates[6] [0]; \ + J[4] = vertex_coordinates[13][0] - vertex_coordinates[12][0]; \ + J[5] = vertex_coordinates[14][0] - vertex_coordinates[12][0]; + +/* Compute Jacobian J for tetrahedron embedded in R^3 */ +#define compute_jacobian_tetrahedron_int_3d(J, vertex_coordinates) \ + J[0] = vertex_coordinates[1] [0] - vertex_coordinates[0] [0]; \ + J[1] = vertex_coordinates[2] [0] - vertex_coordinates[0] [0]; \ + J[2] = vertex_coordinates[3] [0] - vertex_coordinates[0] [0]; \ + J[3] = vertex_coordinates[9] [0] - vertex_coordinates[8] [0]; \ + J[4] = vertex_coordinates[10][0] - vertex_coordinates[8] [0]; \ + J[5] = vertex_coordinates[11][0] - vertex_coordinates[8] [0]; \ + J[6] = vertex_coordinates[17][0] - vertex_coordinates[16][0]; \ + J[7] = vertex_coordinates[18][0] - vertex_coordinates[16][0]; \ + J[8] = vertex_coordinates[19][0] - vertex_coordinates[16][0]; + +/* Compute Jacobian J for tensor product prism embedded in R^3 */ +#define compute_jacobian_prism_int_3d(J, vertex_coordinates) \ + J[0] = vertex_coordinates[2] [0] - vertex_coordinates[0] [0]; \ + J[1] = vertex_coordinates[4] [0] - vertex_coordinates[0] [0]; \ + J[2] = vertex_coordinates[1] [0] - vertex_coordinates[0] [0]; \ + J[3] = vertex_coordinates[14][0] - vertex_coordinates[12][0]; \ + J[4] = vertex_coordinates[16][0] - vertex_coordinates[12][0]; \ + J[5] = vertex_coordinates[13][0] - vertex_coordinates[12][0]; \ + J[6] = vertex_coordinates[26][0] - vertex_coordinates[24][0]; \ + J[7] = vertex_coordinates[28][0] - vertex_coordinates[24][0]; \ + J[8] = vertex_coordinates[25][0] - vertex_coordinates[24][0]; + +/* --- Computation of Jacobian inverses --- */ + +/* Compute Jacobian inverse K for interval embedded in R^1 */ +#define compute_jacobian_inverse_interval_1d(K, det, J) \ + det = J[0]; \ + K[0] = 1.0 / det; + +/* Compute Jacobian (pseudo)inverse K for interval embedded in R^2 */ +#define compute_jacobian_inverse_interval_2d(K, det, J) \ + do { const double det2 = J[0]*J[0] + J[1]*J[1]; \ + det = sqrt(det2); \ + K[0] = J[0] / det2; \ + K[1] = J[1] / det2; } while (0) + +/* Compute Jacobian (pseudo)inverse K for interval embedded in R^3 */ +#define compute_jacobian_inverse_interval_3d(K, det, J) \ + do { const double det2 = J[0]*J[0] + J[1]*J[1] + J[2]*J[2]; \ + det = sqrt(det2); \ + K[0] = J[0] / det2; \ + K[1] = J[1] / det2; \ + K[2] = J[2] / det2; } while (0) + +/* Compute Jacobian inverse K for triangle embedded in R^2 */ +#define compute_jacobian_inverse_triangle_2d(K, det, J) \ + det = J[0]*J[3] - J[1]*J[2]; \ + K[0] = J[3] / det; \ + K[1] = -J[1] / det; \ + K[2] = -J[2] / det; \ + K[3] = J[0] / det; + +/* Compute Jacobian (pseudo)inverse K for triangle embedded in R^3 */ +#define compute_jacobian_inverse_triangle_3d(K, det, J) \ + do { const double d_0 = J[2]*J[5] - J[4]*J[3]; \ + const double d_1 = J[4]*J[1] - J[0]*J[5]; \ + const double d_2 = J[0]*J[3] - J[2]*J[1]; \ + const double c_0 = J[0]*J[0] + J[2]*J[2] + J[4]*J[4]; \ + const double c_1 = J[1]*J[1] + J[3]*J[3] + J[5]*J[5]; \ + const double c_2 = J[0]*J[1] + J[2]*J[3] + J[4]*J[5]; \ + const double den = c_0*c_1 - c_2*c_2; \ + const double det2 = d_0*d_0 + d_1*d_1 + d_2*d_2; \ + det = sqrt(det2); \ + K[0] = (J[0]*c_1 - J[1]*c_2) / den; \ + K[1] = (J[2]*c_1 - J[3]*c_2) / den; \ + K[2] = (J[4]*c_1 - J[5]*c_2) / den; \ + K[3] = (J[1]*c_0 - J[0]*c_2) / den; \ + K[4] = (J[3]*c_0 - J[2]*c_2) / den; \ + K[5] = (J[5]*c_0 - J[4]*c_2) / den; } while (0) + +/* Compute Jacobian (pseudo)inverse K for quad embedded in R^2 */ +#define compute_jacobian_inverse_quad_2d compute_jacobian_inverse_triangle_2d + +/* Compute Jacobian (pseudo)inverse K for quad embedded in R^3 */ +#define compute_jacobian_inverse_quad_3d compute_jacobian_inverse_triangle_3d + +/* Compute Jacobian inverse K for tetrahedron embedded in R^3 */ +#define compute_jacobian_inverse_tetrahedron_3d(K, det, J) \ + do { const double d_00 = J[4]*J[8] - J[5]*J[7]; \ + const double d_01 = J[5]*J[6] - J[3]*J[8]; \ + const double d_02 = J[3]*J[7] - J[4]*J[6]; \ + const double d_10 = J[2]*J[7] - J[1]*J[8]; \ + const double d_11 = J[0]*J[8] - J[2]*J[6]; \ + const double d_12 = J[1]*J[6] - J[0]*J[7]; \ + const double d_20 = J[1]*J[5] - J[2]*J[4]; \ + const double d_21 = J[2]*J[3] - J[0]*J[5]; \ + const double d_22 = J[0]*J[4] - J[1]*J[3]; \ + det = J[0]*d_00 + J[3]*d_10 + J[6]*d_20; \ + K[0] = d_00 / det; \ + K[1] = d_10 / det; \ + K[2] = d_20 / det; \ + K[3] = d_01 / det; \ + K[4] = d_11 / det; \ + K[5] = d_21 / det; \ + K[6] = d_02 / det; \ + K[7] = d_12 / det; \ + K[8] = d_22 / det; } while(0) + +/* Compute Jacobian inverse K for tensor product prism embedded in R^3 - identical to t et */ +#define compute_jacobian_inverse_prism_3d(K, det, J) \ + do { const double d_00 = J[4]*J[8] - J[5]*J[7]; \ + const double d_01 = J[5]*J[6] - J[3]*J[8]; \ + const double d_02 = J[3]*J[7] - J[4]*J[6]; \ + const double d_10 = J[2]*J[7] - J[1]*J[8]; \ + const double d_11 = J[0]*J[8] - J[2]*J[6]; \ + const double d_12 = J[1]*J[6] - J[0]*J[7]; \ + const double d_20 = J[1]*J[5] - J[2]*J[4]; \ + const double d_21 = J[2]*J[3] - J[0]*J[5]; \ + const double d_22 = J[0]*J[4] - J[1]*J[3]; \ + det = J[0]*d_00 + J[3]*d_10 + J[6]*d_20; \ + K[0] = d_00 / det; \ + K[1] = d_10 / det; \ + K[2] = d_20 / det; \ + K[3] = d_01 / det; \ + K[4] = d_11 / det; \ + K[5] = d_21 / det; \ + K[6] = d_02 / det; \ + K[7] = d_12 / det; \ + K[8] = d_22 / det; } while (0) + +/* --- Compute facet edge lengths --- */ + +#define compute_facet_edge_length_tetrahedron_3d(facet, vertex_coordinates) \ + const unsigned int tetrahedron_facet_edge_vertices[4][3][2] = { \ + {{2, 3}, {1, 3}, {1, 2}}, \ + {{2, 3}, {0, 3}, {0, 2}}, \ + {{1, 3}, {0, 3}, {0, 1}}, \ + {{1, 2}, {0, 2}, {0, 1}}, \ + }; \ + double edge_lengths_sqr[3]; \ + for (unsigned int edge = 0; edge < 3; ++edge) \ + { \ + const unsigned int vertex0 = tetrahedron_facet_edge_vertices[facet][edge][0]; \ + const unsigned int vertex1 = tetrahedron_facet_edge_vertices[facet][edge][1]; \ + edge_lengths_sqr[edge] = (vertex_coordinates[vertex1 + 0][0] - vertex_coordinates[vertex0 + 0][0])*(vertex_coordinates[vertex1 + 0][0] - vertex_coordinates[vertex0 + 0][0]) \ + + (vertex_coordinates[vertex1 + 4][0] - vertex_coordinates[vertex0 + 4][0])*(vertex_coordinates[vertex1 + 4][0] - vertex_coordinates[vertex0 + 4][0]) \ + + (vertex_coordinates[vertex1 + 8][0] - vertex_coordinates[vertex0 + 8][0])*(vertex_coordinates[vertex1 + 8][0] - vertex_coordinates[vertex0 + 8][0]); \ + } + +/* Compute min edge length in facet of tetrahedron embedded in R^3 */ +#define compute_min_facet_edge_length_tetrahedron_3d(min_edge_length, facet, vertex_coordinates) \ + compute_facet_edge_length_tetrahedron_3d(facet, vertex_coordinates); \ + min_edge_length = sqrt(fmin(fmin(edge_lengths_sqr[1], edge_lengths_sqr[1]), edge_lengths_sqr[2])); + +/* Compute max edge length in facet of tetrahedron embedded in R^3 */ +/* + * FIXME: we can't call compute_facet_edge_length_tetrahedron_3d again, so we + * rely on the fact that max is always computed after min + */ +#define compute_max_facet_edge_length_tetrahedron_3d(max_edge_length, facet, vertex_coordinates) \ + max_edge_length = sqrt(fmax(fmax(edge_lengths_sqr[1], edge_lengths_sqr[1]), edge_lengths_sqr[2])); diff --git a/firedrake/solving.py b/firedrake/solving.py index 4072e322f2..5fc69dc56f 100644 --- a/firedrake/solving.py +++ b/firedrake/solving.py @@ -28,11 +28,12 @@ import ufl from ufl_expr import derivative -from pyop2 import op2, ffc_interface +from pyop2 import op2 from pyop2.exceptions import MapValueError from pyop2.logger import progress, INFO import core_types import types +from ffc_interface import compile_form from assemble_expressions import assemble_expression from petsc import PETSc @@ -380,7 +381,7 @@ def _assemble(f, tensor=None, bcs=None): """ - kernels = ffc_interface.compile_form(f, "form") + kernels = compile_form(f, "form") fd = f.form_data() diff --git a/tests/test_ffc_interface.py b/tests/test_ffc_interface.py new file mode 100644 index 0000000000..4bc17c51a0 --- /dev/null +++ b/tests/test_ffc_interface.py @@ -0,0 +1,106 @@ +import pytest +from firedrake import * +import os + + +@pytest.fixture(scope='module') +def fs(): + mesh = UnitSquareMesh(1, 1) + return FunctionSpace(mesh, 'CG', 1) + + +@pytest.fixture +def mass(fs): + u = TestFunction(fs) + v = TrialFunction(fs) + return u * v * dx + + +@pytest.fixture +def laplace(fs): + u = TestFunction(fs) + v = TrialFunction(fs) + return inner(grad(u), grad(v)) * dx + + +@pytest.fixture +def rhs(fs): + v = TrialFunction(fs) + g = Function(fs) + return g * v * ds + + +@pytest.fixture +def rhs2(fs): + v = TrialFunction(fs) + f = Function(fs) + g = Function(fs) + return f * v * dx + g * v * ds + + +@pytest.fixture +def cache_key(mass): + return ffc_interface.FFCKernel(mass, 'mass').cache_key + + +@pytest.mark.xfail("not hasattr(ffc_interface.constants, 'PYOP2_VERSION')") +class TestFFCCache: + + """FFC code generation cache tests.""" + + def test_ffc_cache_dir_exists(self): + """Importing ffc_interface should create FFC Kernel cache dir.""" + assert os.path.exists(ffc_interface.FFCKernel._cachedir) + + def test_ffc_cache_persist_on_disk(self, cache_key): + """FFCKernel should be persisted on disk.""" + assert os.path.exists( + os.path.join(ffc_interface.FFCKernel._cachedir, cache_key)) + + def test_ffc_cache_read_from_disk(self, cache_key): + """Loading an FFCKernel from disk should yield the right object.""" + assert ffc_interface.FFCKernel._read_from_disk( + cache_key).cache_key == cache_key + + def test_ffc_compute_form_data(self, mass): + """Compiling a form attaches form data.""" + ffc_interface.compile_form(mass, 'mass') + + assert mass.form_data() + + def test_ffc_same_form(self, mass): + """Compiling the same form twice should load kernels from cache.""" + k1 = ffc_interface.compile_form(mass, 'mass') + k2 = ffc_interface.compile_form(mass, 'mass') + + assert k1 is k2 + + def test_ffc_different_forms(self, mass, laplace): + """Compiling different forms should not load kernels from cache.""" + k1 = ffc_interface.compile_form(mass, 'mass') + k2 = ffc_interface.compile_form(laplace, 'mass') + + assert k1 is not k2 + + def test_ffc_different_names(self, mass): + """Compiling different forms should not load kernels from cache.""" + k1 = ffc_interface.compile_form(mass, 'mass') + k2 = ffc_interface.compile_form(mass, 'laplace') + + assert k1 is not k2 + + def test_ffc_cell_kernel(self, mass): + k = ffc_interface.compile_form(mass, 'mass') + assert 'cell_integral' in k[0].code and len(k) == 1 + + def test_ffc_exterior_facet_kernel(self, rhs): + k = ffc_interface.compile_form(rhs, 'rhs') + assert 'exterior_facet_integral' in k[0].code and len(k) == 1 + + def test_ffc_cell_exterior_facet_kernel(self, rhs2): + k = ffc_interface.compile_form(rhs2, 'rhs2') + assert 'cell_integral' in k[ + 0].code and 'exterior_facet_integral' in k[1].code and len(k) == 2 + +if __name__ == '__main__': + pytest.main(os.path.abspath(__file__)) From 91beb3d1cb346c62be7a2ca3db578d107e65d033 Mon Sep 17 00:00:00 2001 From: Florian Rathgeber Date: Fri, 31 Jan 2014 10:17:29 +0000 Subject: [PATCH 02/27] Add compatible FFC version --- firedrake/ffc_interface.py | 2 +- firedrake/version.py | 2 ++ 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/firedrake/ffc_interface.py b/firedrake/ffc_interface.py index 6fbfabd61c..3e0a9962bc 100644 --- a/firedrake/ffc_interface.py +++ b/firedrake/ffc_interface.py @@ -28,7 +28,7 @@ def _check_version(): - from pyop2.version import __compatible_ffc_version_info__ as compatible_version, \ + from version import __compatible_ffc_version_info__ as compatible_version, \ __compatible_ffc_version__ as version try: if constants.PYOP2_VERSION_INFO[:2] == compatible_version[:2]: diff --git a/firedrake/version.py b/firedrake/version.py index df98e796d0..a7e848909d 100644 --- a/firedrake/version.py +++ b/firedrake/version.py @@ -1,5 +1,7 @@ __version_info__ = (0, 10, 0) __version__ = '.'.join(map(str, __version_info__)) +__compatible_ffc_version_info__ = (0, 5, 0) +__compatible_ffc_version__ = '.'.join(map(str, __compatible_ffc_version_info__)) def check(): From 6324bd8c73386f40c257e6bc9cacc25af11c15a6 Mon Sep 17 00:00:00 2001 From: Florian Rathgeber Date: Mon, 3 Feb 2014 14:32:38 +0000 Subject: [PATCH 03/27] First stab at a FormSplitter --- firedrake/ffc_interface.py | 121 ++++++++++++++++++++++++++++++++++++- 1 file changed, 119 insertions(+), 2 deletions(-) diff --git a/firedrake/ffc_interface.py b/firedrake/ffc_interface.py index 3e0a9962bc..224ec59d9c 100644 --- a/firedrake/ffc_interface.py +++ b/firedrake/ffc_interface.py @@ -4,9 +4,13 @@ from hashlib import md5 import os import tempfile +import numpy as np + +from ufl import Form, FiniteElement, VectorElement +from ufl.algorithms import as_form, traverse_terminals, ReuseTransformer +from ufl.indexing import FixedIndex, MultiIndex +from ufl_expr import Argument -from ufl import Form -from ufl.algorithms import as_form from ffc import default_parameters, compile_form as ffc_compile_form from ffc import constants @@ -14,6 +18,7 @@ from pyop2.op2 import Kernel from pyop2.mpi import MPI from pyop2.ir.ast_base import PreprocessNode, Root +from pyop2.utils import as_tuple _form_cache = {} @@ -39,6 +44,118 @@ def _check_version(): % (version, getattr(constants, 'PYOP2_VERSION', 'unknown'))) +class FormSplitter(ReuseTransformer): + """Split a form into a subtree for each component of the mixed space it is + built on. This is a no-op on forms over non-mixed spaces.""" + + def split(self, form): + """Split the given form.""" + fd = form.compute_form_data() + # If there is no mixed element involved, return a form per integral + if all(isinstance(e, (FiniteElement, VectorElement)) for e in fd.unique_sub_elements): + return [[Form([i])] for i in form.integrals()] + # Otherwise visit each integrand and obtain the tuple of sub forms + return [[f * i.measure() for f in as_tuple(self.visit(i.integrand()))] + for i in form.integrals()] + + def sum(self, o, l, r): + """Take the sum of operands on the same block and return a tuple of + partial sums for each block.""" + + def find_idx(e): + """Find the block index of an expression given by the indices of + the function spaces of the arguments (test and trial function).""" + row, col = None, None + for t in traverse_terminals(e): + if isinstance(t, Argument): + if t.count() == -2: # Test function gives the row + row = t.function_space().index + elif t.count() == -1: # Trial function gives the column + col = t.function_space().index + return (row, col) + + as_list = lambda o: list(o) if isinstance(o, (list, tuple)) else [o] + res = [] + # For each (index, argument) tuple in the left operand list, look for + # a tuple with corresponding index in the right operand list. If + # there is one, append the sum of the arguments with that index to the + # results list, otherwise just the tuple from the left operand list + l = as_list(l) + r = as_list(r) + idx_r = [find_idx(i) for i in r] + # Go over all the operands in the left operand list + for a, i in zip(l, [find_idx(i) for i in l]): + # If there is any operand in the right operand list on the same + # block, take their sum + try: + j = idx_r.index(i) + idx_r.pop(j) + res.append(o.reconstruct(a, r.pop(j))) + # Otherwise just append the operand from the left operand list + except ValueError: + res.append(a) + # All remaining tuples in the right operand list had no matches, so we + # append them to the results list + return tuple(res + r) + + def _binop(self, o, l, r): + if isinstance(l, tuple) and isinstance(r, tuple): + return tuple(o.reconstruct(op1, op2) for op1, op2 in zip(l, r)) + else: + return o.reconstruct(l, r) + + def inner(self, o, l, r): + """Reconstruct an inner product on each of the component spaces.""" + return self._binop(o, l, r) + + def product(self, o, l, r): + """Reconstruct a product on each of the component spaces.""" + return self._binop(o, l, r) + + def dot(self, o, l, r): + """Reconstruct a dot product on each of the component spaces.""" + return self._binop(o, l, r) + + def indexed(self, o, arg, idx): + """Apply fixed indices where they point on a scalar subspace. + Reconstruct fixed indices on a component vector and any other index.""" + if isinstance(idx._indices[0], FixedIndex): + # Find the element to which the FixedIndex points. We might deal + # with coefficients on vector elements, in which case we need to + # reconstruct the indexed with an adjusted index space. Otherwise + # we can just return the coefficient. + i = idx._indices[0]._value + pos = 0 + for op in arg: + # If the FixedIndex points at a scalar (shapeless) operand, + # return it + if not op.shape() and i == pos: + return op + size = np.prod(op.shape() or 1) + # If the FixedIndex points at a component of the current + # operand, reconstruct an Indexed with an adjusted index space + if i < pos + size: + return o.reconstruct(op, MultiIndex(FixedIndex(i - pos), {})) + # Otherwise update the position in the index space + pos += size + raise NotImplementedError("No idea what to in %r with %r" % (o, arg)) + else: + return o.reconstruct(arg, idx) + + def argument(self, o): + """Split an argument into its constituent spaces.""" + if isinstance(o.element(), (FiniteElement, VectorElement)): + return o + return tuple(Argument(fs.ufl_element(), fs, o.count()) + for fs in o.function_space().split()) + + def coefficient(self, o): + """Split a coefficient into its constituent spaces.""" + if isinstance(o.element(), (FiniteElement, VectorElement)): + return o + return o.split() + + class FFCKernel(DiskCached): _cache = {} From 14f41099368d450477a568eebea9d570707c0eb5 Mon Sep 17 00:00:00 2001 From: Florian Rathgeber Date: Tue, 4 Feb 2014 12:20:07 +0000 Subject: [PATCH 04/27] Sum integrands on the same measure before splitting a form --- firedrake/ffc_interface.py | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/firedrake/ffc_interface.py b/firedrake/ffc_interface.py index 224ec59d9c..8937906c9c 100644 --- a/firedrake/ffc_interface.py +++ b/firedrake/ffc_interface.py @@ -2,6 +2,7 @@ generated code in order to make it suitable for passing to the backends.""" from hashlib import md5 +from operator import add import os import tempfile import numpy as np @@ -44,6 +45,12 @@ def _check_version(): % (version, getattr(constants, 'PYOP2_VERSION', 'unknown'))) +def sum_integrands(form): + """Produce a form with the integrands on the same measure summed.""" + return Form([it[0].reconstruct(reduce(add, [i.integrand() for i in it])) + for d, it in form.integral_groups().items()]) + + class FormSplitter(ReuseTransformer): """Split a form into a subtree for each component of the mixed space it is built on. This is a no-op on forms over non-mixed spaces.""" @@ -53,10 +60,10 @@ def split(self, form): fd = form.compute_form_data() # If there is no mixed element involved, return a form per integral if all(isinstance(e, (FiniteElement, VectorElement)) for e in fd.unique_sub_elements): - return [[Form([i])] for i in form.integrals()] + return [[Form([i])] for i in sum_integrands(form).integrals()] # Otherwise visit each integrand and obtain the tuple of sub forms return [[f * i.measure() for f in as_tuple(self.visit(i.integrand()))] - for i in form.integrals()] + for i in sum_integrands(form).integrals()] def sum(self, o, l, r): """Take the sum of operands on the same block and return a tuple of From 11211f08faa85b9c02d000e93fbfd4559b7afe55 Mon Sep 17 00:00:00 2001 From: Florian Rathgeber Date: Tue, 4 Feb 2014 13:52:13 +0000 Subject: [PATCH 05/27] Reconstruct Indexed and IndexSum for non fixed indices --- firedrake/ffc_interface.py | 16 ++++++++++++++-- 1 file changed, 14 insertions(+), 2 deletions(-) diff --git a/firedrake/ffc_interface.py b/firedrake/ffc_interface.py index 8937906c9c..90c2c7cb2b 100644 --- a/firedrake/ffc_interface.py +++ b/firedrake/ffc_interface.py @@ -123,6 +123,19 @@ def dot(self, o, l, r): """Reconstruct a dot product on each of the component spaces.""" return self._binop(o, l, r) + def _index(self, o, arg, idx): + """Reconstruct an index if the rank matches, otherwise yield the + argument. If the argument is a tuple, go over each entry.""" + build = lambda a: o.reconstruct(a, idx) if a.rank() == len(idx.free_indices()) else a + if isinstance(arg, tuple): + return tuple(build(a) for a in arg) + else: + return build(arg) + + def index_sum(self, o, arg, idx): + """Reconstruct an index sum on each of the component spaces.""" + return self._index(o, arg, idx) + def indexed(self, o, arg, idx): """Apply fixed indices where they point on a scalar subspace. Reconstruct fixed indices on a component vector and any other index.""" @@ -146,8 +159,7 @@ def indexed(self, o, arg, idx): # Otherwise update the position in the index space pos += size raise NotImplementedError("No idea what to in %r with %r" % (o, arg)) - else: - return o.reconstruct(arg, idx) + return self._index(o, arg, idx) def argument(self, o): """Split an argument into its constituent spaces.""" From 77ea55649d1e05b33ba939c2188b563df2c6b968 Mon Sep 17 00:00:00 2001 From: Florian Rathgeber Date: Mon, 10 Feb 2014 14:20:08 +0000 Subject: [PATCH 06/27] Check for MixedFunctionSpace when splitting Argument/Coefficient --- firedrake/ffc_interface.py | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/firedrake/ffc_interface.py b/firedrake/ffc_interface.py index 90c2c7cb2b..ee8be83d67 100644 --- a/firedrake/ffc_interface.py +++ b/firedrake/ffc_interface.py @@ -21,6 +21,8 @@ from pyop2.ir.ast_base import PreprocessNode, Root from pyop2.utils import as_tuple +import types + _form_cache = {} ffc_parameters = default_parameters() @@ -163,16 +165,16 @@ def indexed(self, o, arg, idx): def argument(self, o): """Split an argument into its constituent spaces.""" - if isinstance(o.element(), (FiniteElement, VectorElement)): - return o - return tuple(Argument(fs.ufl_element(), fs, o.count()) - for fs in o.function_space().split()) + if isinstance(o.function_space(), types.MixedFunctionSpace): + return tuple(Argument(fs.ufl_element(), fs, o.count()) + for fs in o.function_space().split()) + return o def coefficient(self, o): """Split a coefficient into its constituent spaces.""" - if isinstance(o.element(), (FiniteElement, VectorElement)): - return o - return o.split() + if isinstance(o.function_space(), types.MixedFunctionSpace): + return o.split() + return o class FFCKernel(DiskCached): From ccb01557a63148c5a0e5e1dc96df2c403a0f0593 Mon Sep 17 00:00:00 2001 From: Florian Rathgeber Date: Mon, 10 Feb 2014 14:21:28 +0000 Subject: [PATCH 07/27] Don't return a tuple if a Sum only has one block --- firedrake/ffc_interface.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/firedrake/ffc_interface.py b/firedrake/ffc_interface.py index ee8be83d67..9ae0706e17 100644 --- a/firedrake/ffc_interface.py +++ b/firedrake/ffc_interface.py @@ -105,7 +105,7 @@ def find_idx(e): res.append(a) # All remaining tuples in the right operand list had no matches, so we # append them to the results list - return tuple(res + r) + return tuple(res + r) if len(res + r) > 1 else (res + r)[0] def _binop(self, o, l, r): if isinstance(l, tuple) and isinstance(r, tuple): From 4f077f64625078b60a138102a3463461dda73d35 Mon Sep 17 00:00:00 2001 From: Florian Rathgeber Date: Mon, 10 Feb 2014 14:25:46 +0000 Subject: [PATCH 08/27] Fix condition for reconstructing index_sum --- firedrake/ffc_interface.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/firedrake/ffc_interface.py b/firedrake/ffc_interface.py index 9ae0706e17..6b636a3bfb 100644 --- a/firedrake/ffc_interface.py +++ b/firedrake/ffc_interface.py @@ -136,7 +136,11 @@ def _index(self, o, arg, idx): def index_sum(self, o, arg, idx): """Reconstruct an index sum on each of the component spaces.""" - return self._index(o, arg, idx) + build = lambda a: o.reconstruct(a, idx) if len(a.free_indices()) == len(idx.free_indices()) else a + if isinstance(arg, tuple): + return tuple(build(a) for a in arg) + else: + return build(arg) def indexed(self, o, arg, idx): """Apply fixed indices where they point on a scalar subspace. From 6fdea5cc03e6317582eace0b8349628093c57cb1 Mon Sep 17 00:00:00 2001 From: Florian Rathgeber Date: Tue, 11 Feb 2014 13:15:07 +0000 Subject: [PATCH 09/27] Split forms when building FFCKernel Does not work with disk caching for now since Coefficients can't be serialised. --- firedrake/ffc_interface.py | 53 +++++++++++++++++++++++++------------- 1 file changed, 35 insertions(+), 18 deletions(-) diff --git a/firedrake/ffc_interface.py b/firedrake/ffc_interface.py index 6b636a3bfb..7de4932b9b 100644 --- a/firedrake/ffc_interface.py +++ b/firedrake/ffc_interface.py @@ -15,7 +15,7 @@ from ffc import default_parameters, compile_form as ffc_compile_form from ffc import constants -from pyop2.caching import DiskCached +from pyop2.caching import Cached from pyop2.op2 import Kernel from pyop2.mpi import MPI from pyop2.ir.ast_base import PreprocessNode, Root @@ -181,7 +181,7 @@ def coefficient(self, o): return o -class FFCKernel(DiskCached): +class FFCKernel(Cached): _cache = {} _cachedir = os.path.join(tempfile.gettempdir(), @@ -194,34 +194,51 @@ def _cache_key(cls, form, name): _firedrake_geometry_md5 + constants.FFC_VERSION + constants.PYOP2_VERSION).hexdigest() - def __init__(self, form, name): + def __init__(self, original_form, name): if self._initialized: return incl = PreprocessNode('#include "firedrake_geometry.h"\n') inc = [os.path.dirname(__file__)] - forms = ffc_compile_form(form, prefix=name, parameters=ffc_parameters) - fdict = dict((f.name, f) for f in forms) kernels = [] - for ida in form.form_data().preprocessed_form.integrals(): - fname = '%s_%s_integral_0_%s' % (name, ida.domain_type(), ida.domain_id()) - # Set optimization options - opts = {} if ida.domain_type() not in ['cell'] else \ - {'licm': False, - 'tile': None, - 'vect': None, - 'ap': False, - 'split': None} - kernels.append(Kernel(Root([incl, fdict[fname]]), fname, opts, inc)) + # Note that split forms are batched by integral i.e. they will only + # ever contain a single integral. We therefore always return the first + # element of any lists that contain different integrals. + for forms in FormSplitter().split(original_form): + for i, form in enumerate(forms): + tree, = ffc_compile_form(form, prefix=name + str(i), + parameters=ffc_parameters) + + fd = form.form_data() + ida = fd.integral_data[0] + # Set optimization options + opts = {} if ida.domain_type not in ['cell'] else \ + {'licm': False, + 'tile': None, + 'vect': None, + 'ap': False, + 'split': None} + + fname = '%s%d_%s_integral_0_%s' % (name, i, ida.domain_type, + ida.domain_id) + + if len(forms) == 1: + idx = (0, 0) + else: + t = tuple(a.function_space().index or 0 + for a in fd.original_arguments) or (i, 0) + idx = t if len(t) == 2 else t + (0,) * (2 - len(t)) + kernels.append((idx, ida.integrals[0].measure(), + fd.original_coefficients, + Kernel(Root([incl, tree]), fname, opts, inc))) self.kernels = tuple(kernels) - self._initialized = True def compile_form(form, name): - """Compile a form using FFC and return a tuple of - :class:`Kernels `.""" + """Compile a form using FFC and return a tuple of tuples of + (index, domain type, coefficients, :class:`Kernels `).""" # Check that we get a Form if not isinstance(form, Form): From a3786ee198d9407acd22099aced6cfbe6c3d0312 Mon Sep 17 00:00:00 2001 From: Florian Rathgeber Date: Tue, 11 Feb 2014 17:12:58 +0000 Subject: [PATCH 10/27] Use split forms when assembling Adapt to changes in compile_form, which now returns a four-tuple of index, domain type, coefficients and kernel. We now need to pull apart the tensor and test and trial functions if it is a mixed type and extract the component for the block the kernel contributes to and filter the boundary conditions to only include those defined on the current block. --- firedrake/solving.py | 183 +++++++++++++++++++++++-------------------- 1 file changed, 96 insertions(+), 87 deletions(-) diff --git a/firedrake/solving.py b/firedrake/solving.py index 5fc69dc56f..6e9b228b16 100644 --- a/firedrake/solving.py +++ b/firedrake/solving.py @@ -525,18 +525,38 @@ def thunk(bcs): bottom = any(bc.sub_domain == "bottom" for bc in bcs) top = any(bc.sub_domain == "top" for bc in bcs) extruded_bcs = (bottom, top) - for kernel, integral in zip(kernels, integrals): - domain_type = integral.measure().domain_type() - if domain_type == 'cell': + for (i, j), measure, coefficients, kernel in kernels: + # Extract block from tensor and test/trial spaces + # FIXME Ugly variable renaming required because functions are not + # lexical closures in Python and we're writing to these variables + if is_mat and tensor.sparsity.shape > (1, 1): + t = tensor[i, j] + ts = test.function_space()[i] + tr = trial.function_space()[j] + tsbc = [bc for bc in bcs if bc.function_space().index == i] + trbc = [bc for bc in bcs if bc.function_space().index == j] + elif is_mat: + t = tensor + ts, tr = test, trial + tsbc, trbc = bcs, bcs + elif is_vec and len(tensor) > 1: + t = tensor[i] + ts = test.function_space()[i] + elif is_vec: + t = tensor + ts = test + else: + t = tensor + if measure.domain_type() == 'cell': if is_mat: - tensor_arg = tensor(op2.INC, (test.cell_node_map(bcs)[op2.i[0]], - trial.cell_node_map(bcs)[op2.i[1]]), - flatten=has_vec_fs(test)) + tensor_arg = t(op2.INC, (ts.cell_node_map(tsbc)[op2.i[0]], + tr.cell_node_map(trbc)[op2.i[1]]), + flatten=has_vec_fs(test)) elif is_vec: - tensor_arg = tensor(op2.INC, test.cell_node_map()[op2.i[0]], - flatten=has_vec_fs(test)) + tensor_arg = t(op2.INC, ts.cell_node_map()[op2.i[0]], + flatten=has_vec_fs(test)) else: - tensor_arg = tensor(op2.INC) + tensor_arg = t(op2.INC) itspace = m.cell_set itspace._extruded_bcs = extruded_bcs @@ -546,7 +566,7 @@ def thunk(bcs): if needs_orientations: args.append(coords.function_space().mesh()._cell_orientations(op2.READ)) - for c in fd.original_coefficients: + for c in coefficients: args.append(c.dat(op2.READ, c.cell_node_map(), flatten=has_vec_fs(c))) @@ -554,27 +574,26 @@ def thunk(bcs): op2.par_loop(*args) except MapValueError: raise RuntimeError("Integral measure does not match measure of all coefficients/arguments") - elif domain_type == 'exterior_facet': + elif measure.domain_type() == 'exterior_facet': if op2.MPI.parallel: raise \ NotImplementedError( "No support for facet integrals under MPI yet") if is_mat: - tensor_arg = tensor(op2.INC, - (test.exterior_facet_node_map(bcs)[op2.i[0]], - trial.exterior_facet_node_map(bcs)[op2.i[1]]), - flatten=has_vec_fs(test)) + tensor_arg = t(op2.INC, (ts.exterior_facet_node_map(tsbc)[op2.i[0]], + tr.exterior_facet_node_map(trbc)[op2.i[1]]), + flatten=has_vec_fs(test)) elif is_vec: - tensor_arg = tensor(op2.INC, - test.exterior_facet_node_map()[op2.i[0]], - flatten=has_vec_fs(test)) + tensor_arg = t(op2.INC, + ts.exterior_facet_node_map()[op2.i[0]], + flatten=has_vec_fs(test)) else: - tensor_arg = tensor(op2.INC) - args = [kernel, m.exterior_facets.measure_set(integral.measure()), tensor_arg, + tensor_arg = t(op2.INC) + args = [kernel, m.exterior_facets.measure_set(measure), tensor_arg, coords.dat(op2.READ, coords.exterior_facet_node_map(), flatten=True)] - for c in fd.original_coefficients: + for c in coefficients: args.append(c.dat(op2.READ, c.exterior_facet_node_map(), flatten=has_vec_fs(c))) args.append(m.exterior_facets.local_facet_dat(op2.READ)) @@ -583,23 +602,21 @@ def thunk(bcs): except MapValueError: raise RuntimeError("Integral measure does not match measure of all coefficients/arguments") - elif domain_type in ['exterior_facet_top', 'exterior_facet_bottom']: + elif measure.domain_type() in ['exterior_facet_top', 'exterior_facet_bottom']: if op2.MPI.parallel: raise \ NotImplementedError( "No support for facet integrals under MPI yet") if is_mat: - tensor_arg = tensor(op2.INC, - (test.cell_node_map(bcs)[op2.i[0]], - trial.cell_node_map(bcs)[op2.i[1]]), - flatten=has_vec_fs(test)) + tensor_arg = t(op2.INC, (ts.cell_node_map(tsbc)[op2.i[0]], + tr.cell_node_map(trbc)[op2.i[1]]), + flatten=has_vec_fs(test)) elif is_vec: - tensor_arg = tensor(op2.INC, - test.cell_node_map()[op2.i[0]], - flatten=has_vec_fs(test)) + tensor_arg = t(op2.INC, ts.cell_node_map()[op2.i[0]], + flatten=has_vec_fs(test)) else: - tensor_arg = tensor(op2.INC) + tensor_arg = t(op2.INC) #In the case of extruded meshes with horizontal facet integrals, two #parallel loops will (potentially) get created and called based on the @@ -614,7 +631,7 @@ def thunk(bcs): args = [kernel, set, tensor_arg, coords.dat(op2.READ, coords.cell_node_map(), flatten=True)] - for c in fd.original_coefficients: + for c in coefficients: args.append(c.dat(op2.READ, c.cell_node_map(), flatten=has_vec_fs(c))) try: @@ -622,28 +639,27 @@ def thunk(bcs): except MapValueError: raise RuntimeError("Integral measure does not match measure of all coefficients/arguments") - elif domain_type == 'exterior_facet_vert': + elif measure.domain_type() == 'exterior_facet_vert': if op2.MPI.parallel: raise \ NotImplementedError( "No support for facet integrals under MPI yet") if is_mat: - tensor_arg = tensor(op2.INC, - (test.exterior_facet_node_map(bcs)[op2.i[0]], - trial.exterior_facet_node_map(bcs)[op2.i[1]]), - flatten=has_vec_fs(test)) + tensor_arg = t(op2.INC, + (ts.exterior_facet_node_map(tsbc)[op2.i[0]], + tr.exterior_facet_node_map(trbc)[op2.i[1]]), + flatten=has_vec_fs(test)) elif is_vec: - tensor_arg = tensor(op2.INC, - test.exterior_facet_node_map()[op2.i[0]], - flatten=has_vec_fs(test)) + tensor_arg = t(op2.INC, ts.exterior_facet_node_map()[op2.i[0]], + flatten=has_vec_fs(test)) else: - tensor_arg = tensor(op2.INC) + tensor_arg = t(op2.INC) args = [kernel, m.exterior_facets.measure_set(integral.measure()), tensor_arg, coords.dat(op2.READ, coords.exterior_facet_node_map(), flatten=True)] - for c in fd.original_coefficients: + for c in coefficients: args.append(c.dat(op2.READ, c.exterior_facet_node_map(), flatten=has_vec_fs(c))) args.append(m.exterior_facets.local_facet_dat(op2.READ)) @@ -652,28 +668,26 @@ def thunk(bcs): except MapValueError: raise RuntimeError("Integral measure does not match measure of all coefficients/arguments") - elif domain_type == 'interior_facet': + elif measure.domain_type() == 'interior_facet': if op2.MPI.parallel: raise \ NotImplementedError( "No support for facet integrals under MPI yet") if is_mat: - tensor_arg = tensor( - op2.INC, (test.interior_facet_node_map(bcs)[op2.i[0]], - trial.interior_facet_node_map(bcs)[ - op2.i[1]]), - flatten=True) + tensor_arg = t(op2.INC, + (ts.interior_facet_node_map(tsbc)[op2.i[0]], + tr.interior_facet_node_map(trbc)[op2.i[1]]), + flatten=True) elif is_vec: - tensor_arg = tensor( - op2.INC, test.interior_facet_node_map()[op2.i[0]], - flatten=True) + tensor_arg = t(op2.INC, ts.interior_facet_node_map()[op2.i[0]], + flatten=True) else: - tensor_arg = tensor(op2.INC) + tensor_arg = t(op2.INC) args = [kernel, m.interior_facets.set, tensor_arg, coords.dat(op2.READ, coords.interior_facet_node_map(), flatten=True)] - for c in fd.original_coefficients: + for c in coefficients: args.append(c.dat(op2.READ, c.interior_facet_node_map(), flatten=True)) args.append(m.interior_facets.local_facet_dat(op2.READ)) @@ -682,28 +696,26 @@ def thunk(bcs): except MapValueError: raise RuntimeError("Integral measure does not match measure of all coefficients/arguments") - elif domain_type == 'interior_facet_horiz': + elif measure.domain_type() == 'interior_facet_horiz': if op2.MPI.parallel: raise \ NotImplementedError( "No support for facet integrals under MPI yet") if is_mat: - tensor_arg = tensor( - op2.INC, (test.cell_node_map(bcs)[op2.i[0]], - trial.cell_node_map(bcs)[op2.i[1]]), - flatten=True) + tensor_arg = t(op2.INC, (ts.cell_node_map(tsbc)[op2.i[0]], + tr.cell_node_map(trbc)[op2.i[1]]), + flatten=True) elif is_vec: - tensor_arg = tensor( - op2.INC, test.cell_node_map()[op2.i[0]], - flatten=True) + tensor_arg = t(op2.INC, ts.cell_node_map()[op2.i[0]], + flatten=True) else: - tensor_arg = tensor(op2.INC) + tensor_arg = t(op2.INC) args = [kernel, m.interior_facets.measure_set(integral.measure()), tensor_arg, coords.dat(op2.READ, coords.cell_node_map(), flatten=True)] - for c in fd.original_coefficients: + for c in coefficients: args.append(c.dat(op2.READ, c.cell_node_map(), flatten=has_vec_fs(c))) try: @@ -711,28 +723,27 @@ def thunk(bcs): except MapValueError: raise RuntimeError("Integral measure does not match measure of all coefficients/arguments") - elif domain_type == 'interior_facet_vert': + elif measure.domain_type() == 'interior_facet_vert': if op2.MPI.parallel: raise \ NotImplementedError( "No support for facet integrals under MPI yet") if is_mat: - tensor_arg = tensor( - op2.INC, (test.interior_facet_node_map(bcs)[op2.i[0]], - trial.interior_facet_node_map(bcs)[ - op2.i[1]]), - flatten=True) + tensor_arg = t(op2.INC, + (ts.interior_facet_node_map(tsbc)[op2.i[0]], + tr.interior_facet_node_map(trbc)[op2.i[1]]), + flatten=True) elif is_vec: - tensor_arg = tensor( - op2.INC, test.interior_facet_node_map()[op2.i[0]], - flatten=True) + tensor_arg = t(op2.INC, + ts.interior_facet_node_map()[op2.i[0]], + flatten=True) else: - tensor_arg = tensor(op2.INC) + tensor_arg = t(op2.INC) args = [kernel, m.interior_facets.set, tensor_arg, coords.dat(op2.READ, coords.interior_facet_node_map(), flatten=True)] - for c in fd.original_coefficients: + for c in coefficients: args.append(c.dat(op2.READ, c.interior_facet_node_map(), flatten=has_vec_fs(c))) args.append(m.interior_facets.local_facet_dat(op2.READ)) @@ -742,20 +753,18 @@ def thunk(bcs): raise RuntimeError("Integral measure does not match measure of all coefficients/arguments") else: - raise RuntimeError('Unknown domain type "%s"' % domain_type) - - if bcs is not None and is_mat: - for bc in bcs: - fs = bc.function_space() - if isinstance(fs, types.MixedFunctionSpace): - raise RuntimeError("""Cannot apply boundary conditions to full mixed space. Did you forget to index it?""") - # Set diagonal entries on bc nodes to 1. - if fs.index is None: - # Non-mixed case - tensor.inc_local_diagonal_entries(bc.nodes) - else: - # Mixed case with indexed FS, zero appropriate block - tensor[fs.index, fs.index].inc_local_diagonal_entries(bc.nodes) + raise RuntimeError('Unknown domain type "%s"' % measure.domain_type()) + + if bcs is not None and is_mat: + for bc in bcs: + fs = bc.function_space() + if isinstance(fs, types.MixedFunctionSpace): + raise RuntimeError("""Cannot apply boundary conditions to full mixed space. Did you forget to index it?""") + # Set diagonal entries on bc nodes to 1 if the current + # block is on the matrix diagonal and its index matches the + # index of the function space the bc is defined on. + if i == j and (fs.index is None or fs.index == i): + t.inc_local_diagonal_entries(bc.nodes) return result() From d1fb01d868f8149c3958821fd11f637399351954 Mon Sep 17 00:00:00 2001 From: Florian Rathgeber Date: Tue, 11 Feb 2014 17:14:28 +0000 Subject: [PATCH 11/27] Split the form in compile_form FFCKernel goes back to handling only a single kernel on a single integral. --- firedrake/ffc_interface.py | 64 +++++++++++++++++-------------------- tests/test_ffc_interface.py | 26 +++++++-------- 2 files changed, 43 insertions(+), 47 deletions(-) diff --git a/firedrake/ffc_interface.py b/firedrake/ffc_interface.py index 7de4932b9b..33953800e8 100644 --- a/firedrake/ffc_interface.py +++ b/firedrake/ffc_interface.py @@ -15,7 +15,7 @@ from ffc import default_parameters, compile_form as ffc_compile_form from ffc import constants -from pyop2.caching import Cached +from pyop2.caching import DiskCached from pyop2.op2 import Kernel from pyop2.mpi import MPI from pyop2.ir.ast_base import PreprocessNode, Root @@ -181,7 +181,7 @@ def coefficient(self, o): return o -class FFCKernel(Cached): +class FFCKernel(DiskCached): _cache = {} _cachedir = os.path.join(tempfile.gettempdir(), @@ -194,44 +194,25 @@ def _cache_key(cls, form, name): _firedrake_geometry_md5 + constants.FFC_VERSION + constants.PYOP2_VERSION).hexdigest() - def __init__(self, original_form, name): + def __init__(self, form, name): if self._initialized: return incl = PreprocessNode('#include "firedrake_geometry.h"\n') inc = [os.path.dirname(__file__)] + ffc_tree = ffc_compile_form(form, prefix=name, parameters=ffc_parameters) kernels = [] - # Note that split forms are batched by integral i.e. they will only - # ever contain a single integral. We therefore always return the first - # element of any lists that contain different integrals. - for forms in FormSplitter().split(original_form): - for i, form in enumerate(forms): - tree, = ffc_compile_form(form, prefix=name + str(i), - parameters=ffc_parameters) - - fd = form.form_data() - ida = fd.integral_data[0] - # Set optimization options - opts = {} if ida.domain_type not in ['cell'] else \ - {'licm': False, - 'tile': None, - 'vect': None, - 'ap': False, - 'split': None} - - fname = '%s%d_%s_integral_0_%s' % (name, i, ida.domain_type, - ida.domain_id) - - if len(forms) == 1: - idx = (0, 0) - else: - t = tuple(a.function_space().index or 0 - for a in fd.original_arguments) or (i, 0) - idx = t if len(t) == 2 else t + (0,) * (2 - len(t)) - kernels.append((idx, ida.integrals[0].measure(), - fd.original_coefficients, - Kernel(Root([incl, tree]), fname, opts, inc))) + for ida, kernel in zip(form.form_data().integral_data, ffc_tree): + # Set optimization options + opts = {} if ida.domain_type not in ['cell'] else \ + {'licm': False, + 'tile': None, + 'vect': None, + 'ap': False, + 'split': None} + kernels.append(Kernel(Root([incl, kernel]), '%s_%s_integral_0_%s' % + (name, ida.domain_type, ida.domain_id), opts, inc)) self.kernels = tuple(kernels) self._initialized = True @@ -244,7 +225,22 @@ def compile_form(form, name): if not isinstance(form, Form): form = as_form(form) - return FFCKernel(form, name).kernels + kernels = [] + for forms in FormSplitter().split(form): + for i, form in enumerate(forms): + kernel, = FFCKernel(form, name + str(i)).kernels + + fd = form.form_data() + ida = fd.integral_data[0] + if len(forms) == 1 and fd.rank == 0: + idx = (0, 0) + else: + t = tuple(a.function_space().index or 0 + for a in fd.original_arguments) or (i, 0) + idx = t if len(t) == 2 else t + (0,) * (2 - len(t)) + kernels.append((idx, ida.integrals[0].measure(), + fd.original_coefficients, kernel)) + return kernels def clear_cache(): diff --git a/tests/test_ffc_interface.py b/tests/test_ffc_interface.py index 4bc17c51a0..60b46f66d6 100644 --- a/tests/test_ffc_interface.py +++ b/tests/test_ffc_interface.py @@ -70,37 +70,37 @@ def test_ffc_compute_form_data(self, mass): def test_ffc_same_form(self, mass): """Compiling the same form twice should load kernels from cache.""" - k1 = ffc_interface.compile_form(mass, 'mass') - k2 = ffc_interface.compile_form(mass, 'mass') + k1, = ffc_interface.compile_form(mass, 'mass') + k2, = ffc_interface.compile_form(mass, 'mass') - assert k1 is k2 + assert k1[-1] is k2[-1] def test_ffc_different_forms(self, mass, laplace): """Compiling different forms should not load kernels from cache.""" - k1 = ffc_interface.compile_form(mass, 'mass') - k2 = ffc_interface.compile_form(laplace, 'mass') + k1, = ffc_interface.compile_form(mass, 'mass') + k2, = ffc_interface.compile_form(laplace, 'mass') - assert k1 is not k2 + assert k1[-1] is not k2[-1] def test_ffc_different_names(self, mass): """Compiling different forms should not load kernels from cache.""" - k1 = ffc_interface.compile_form(mass, 'mass') - k2 = ffc_interface.compile_form(mass, 'laplace') + k1, = ffc_interface.compile_form(mass, 'mass') + k2, = ffc_interface.compile_form(mass, 'laplace') - assert k1 is not k2 + assert k1[-1] is not k2[-1] def test_ffc_cell_kernel(self, mass): k = ffc_interface.compile_form(mass, 'mass') - assert 'cell_integral' in k[0].code and len(k) == 1 + assert 'cell_integral' in k[0][-1].code and len(k) == 1 def test_ffc_exterior_facet_kernel(self, rhs): k = ffc_interface.compile_form(rhs, 'rhs') - assert 'exterior_facet_integral' in k[0].code and len(k) == 1 + assert 'exterior_facet_integral' in k[0][-1].code and len(k) == 1 def test_ffc_cell_exterior_facet_kernel(self, rhs2): k = ffc_interface.compile_form(rhs2, 'rhs2') - assert 'cell_integral' in k[ - 0].code and 'exterior_facet_integral' in k[1].code and len(k) == 2 + assert 'cell_integral' in k[0][-1].code and \ + 'exterior_facet_integral' in k[1][-1].code and len(k) == 2 if __name__ == '__main__': pytest.main(os.path.abspath(__file__)) From daa105b8f295f653e8891361d8ad0045ae832b13 Mon Sep 17 00:00:00 2001 From: Florian Rathgeber Date: Tue, 11 Feb 2014 17:42:51 +0000 Subject: [PATCH 12/27] Handle grad in FormSplitter --- firedrake/ffc_interface.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/firedrake/ffc_interface.py b/firedrake/ffc_interface.py index 33953800e8..87a49831a8 100644 --- a/firedrake/ffc_interface.py +++ b/firedrake/ffc_interface.py @@ -125,6 +125,13 @@ def dot(self, o, l, r): """Reconstruct a dot product on each of the component spaces.""" return self._binop(o, l, r) + def grad(self, o, arg): + """Reconstruct a grad on each of the component spaces.""" + if isinstance(arg, tuple): + return tuple(o.reconstruct(a) for a in arg) + else: + return o.reconstruct(arg) + def _index(self, o, arg, idx): """Reconstruct an index if the rank matches, otherwise yield the argument. If the argument is a tuple, go over each entry.""" From 5b9a7a586a87fab043fe0170527bb7bff51dfce2 Mon Sep 17 00:00:00 2001 From: Florian Rathgeber Date: Thu, 13 Feb 2014 14:38:30 +0000 Subject: [PATCH 13/27] Simplify FormSplitter using splitting by construction Instead of making one pass over the form and getting back a list of subforms on the blocks of the mixed space, make as many passes as there are blocks in the mixed space and only enabling one component of each Argument at a time. When splitting an Argument, create a UFL ListTensor where only the component of the MixedFunctionSpace selected by the current block index is included and all other components are set to Zero. These Zeros are subsequently eliminated as they are used in expressions and not propagated up the form tree. --- firedrake/ffc_interface.py | 164 ++++++++----------------------------- 1 file changed, 35 insertions(+), 129 deletions(-) diff --git a/firedrake/ffc_interface.py b/firedrake/ffc_interface.py index 87a49831a8..bb85a67f3f 100644 --- a/firedrake/ffc_interface.py +++ b/firedrake/ffc_interface.py @@ -5,11 +5,10 @@ from operator import add import os import tempfile -import numpy as np -from ufl import Form, FiniteElement, VectorElement -from ufl.algorithms import as_form, traverse_terminals, ReuseTransformer -from ufl.indexing import FixedIndex, MultiIndex +from ufl import Form, FiniteElement, VectorElement, as_vector +from ufl.algorithms import as_form, ReuseTransformer +from ufl.constantvalue import Zero from ufl_expr import Argument from ffc import default_parameters, compile_form as ffc_compile_form @@ -19,7 +18,6 @@ from pyop2.op2 import Kernel from pyop2.mpi import MPI from pyop2.ir.ast_base import PreprocessNode, Root -from pyop2.utils import as_tuple import types @@ -62,129 +60,44 @@ def split(self, form): fd = form.compute_form_data() # If there is no mixed element involved, return a form per integral if all(isinstance(e, (FiniteElement, VectorElement)) for e in fd.unique_sub_elements): - return [[Form([i])] for i in sum_integrands(form).integrals()] + return [[((0, 0), Form([i]))] for i in sum_integrands(form).integrals()] # Otherwise visit each integrand and obtain the tuple of sub forms - return [[f * i.measure() for f in as_tuple(self.visit(i.integrand()))] - for i in sum_integrands(form).integrals()] - - def sum(self, o, l, r): - """Take the sum of operands on the same block and return a tuple of - partial sums for each block.""" - - def find_idx(e): - """Find the block index of an expression given by the indices of - the function spaces of the arguments (test and trial function).""" - row, col = None, None - for t in traverse_terminals(e): - if isinstance(t, Argument): - if t.count() == -2: # Test function gives the row - row = t.function_space().index - elif t.count() == -1: # Trial function gives the column - col = t.function_space().index - return (row, col) - - as_list = lambda o: list(o) if isinstance(o, (list, tuple)) else [o] - res = [] - # For each (index, argument) tuple in the left operand list, look for - # a tuple with corresponding index in the right operand list. If - # there is one, append the sum of the arguments with that index to the - # results list, otherwise just the tuple from the left operand list - l = as_list(l) - r = as_list(r) - idx_r = [find_idx(i) for i in r] - # Go over all the operands in the left operand list - for a, i in zip(l, [find_idx(i) for i in l]): - # If there is any operand in the right operand list on the same - # block, take their sum - try: - j = idx_r.index(i) - idx_r.pop(j) - res.append(o.reconstruct(a, r.pop(j))) - # Otherwise just append the operand from the left operand list - except ValueError: - res.append(a) - # All remaining tuples in the right operand list had no matches, so we - # append them to the results list - return tuple(res + r) if len(res + r) > 1 else (res + r)[0] - - def _binop(self, o, l, r): - if isinstance(l, tuple) and isinstance(r, tuple): - return tuple(o.reconstruct(op1, op2) for op1, op2 in zip(l, r)) - else: - return o.reconstruct(l, r) - - def inner(self, o, l, r): - """Reconstruct an inner product on each of the component spaces.""" - return self._binop(o, l, r) - - def product(self, o, l, r): - """Reconstruct a product on each of the component spaces.""" - return self._binop(o, l, r) - - def dot(self, o, l, r): - """Reconstruct a dot product on each of the component spaces.""" - return self._binop(o, l, r) - - def grad(self, o, arg): - """Reconstruct a grad on each of the component spaces.""" - if isinstance(arg, tuple): - return tuple(o.reconstruct(a) for a in arg) - else: - return o.reconstruct(arg) - - def _index(self, o, arg, idx): - """Reconstruct an index if the rank matches, otherwise yield the - argument. If the argument is a tuple, go over each entry.""" - build = lambda a: o.reconstruct(a, idx) if a.rank() == len(idx.free_indices()) else a - if isinstance(arg, tuple): - return tuple(build(a) for a in arg) - else: - return build(arg) - - def index_sum(self, o, arg, idx): - """Reconstruct an index sum on each of the component spaces.""" - build = lambda a: o.reconstruct(a, idx) if len(a.free_indices()) == len(idx.free_indices()) else a - if isinstance(arg, tuple): - return tuple(build(a) for a in arg) - else: - return build(arg) - - def indexed(self, o, arg, idx): - """Apply fixed indices where they point on a scalar subspace. - Reconstruct fixed indices on a component vector and any other index.""" - if isinstance(idx._indices[0], FixedIndex): - # Find the element to which the FixedIndex points. We might deal - # with coefficients on vector elements, in which case we need to - # reconstruct the indexed with an adjusted index space. Otherwise - # we can just return the coefficient. - i = idx._indices[0]._value - pos = 0 - for op in arg: - # If the FixedIndex points at a scalar (shapeless) operand, - # return it - if not op.shape() and i == pos: - return op - size = np.prod(op.shape() or 1) - # If the FixedIndex points at a component of the current - # operand, reconstruct an Indexed with an adjusted index space - if i < pos + size: - return o.reconstruct(op, MultiIndex(FixedIndex(i - pos), {})) - # Otherwise update the position in the index space - pos += size - raise NotImplementedError("No idea what to in %r with %r" % (o, arg)) - return self._index(o, arg, idx) + shape = tuple(len(a.function_space()) for a in fd.original_arguments) + forms_list = [] + for it in sum_integrands(form).integrals(): + forms = [] + for i in range(shape[0] if len(shape) > 0 else 1): + for j in range(shape[1] if len(shape) > 1 else 1): + self._idx = {-2: i, -1: j} + integrand = self.visit(it.integrand()) + if not isinstance(integrand, Zero): + forms.append([((i, j), integrand * it.measure())]) + forms_list += forms + return forms_list def argument(self, o): """Split an argument into its constituent spaces.""" if isinstance(o.function_space(), types.MixedFunctionSpace): - return tuple(Argument(fs.ufl_element(), fs, o.count()) - for fs in o.function_space().split()) + args = [] + for i, fs in enumerate(o.function_space().split()): + a = Argument(fs.ufl_element(), fs, o.count()) + if a.shape(): + if self._idx[o.count()] == i: + args += [a[j] for j in range(a.shape()[0])] + else: + args += [Zero() for j in range(a.shape()[0])] + else: + if self._idx[o.count()] == i: + args.append(a) + else: + args.append(Zero()) + return as_vector(args) return o def coefficient(self, o): """Split a coefficient into its constituent spaces.""" if isinstance(o.function_space(), types.MixedFunctionSpace): - return o.split() + return as_vector(list(o.split())) return o @@ -234,18 +147,11 @@ def compile_form(form, name): kernels = [] for forms in FormSplitter().split(form): - for i, form in enumerate(forms): - kernel, = FFCKernel(form, name + str(i)).kernels - - fd = form.form_data() + for (i, j), form in forms: + fd = form.compute_form_data() + kernel, = FFCKernel(form, name + str(i) + str(j)).kernels ida = fd.integral_data[0] - if len(forms) == 1 and fd.rank == 0: - idx = (0, 0) - else: - t = tuple(a.function_space().index or 0 - for a in fd.original_arguments) or (i, 0) - idx = t if len(t) == 2 else t + (0,) * (2 - len(t)) - kernels.append((idx, ida.integrals[0].measure(), + kernels.append(((i, j), ida.integrals[0].measure(), fd.original_coefficients, kernel)) return kernels From 7d5f3e847ee246f85e0a40a20729702cb236cc42 Mon Sep 17 00:00:00 2001 From: Florian Rathgeber Date: Fri, 14 Feb 2014 14:02:10 +0000 Subject: [PATCH 14/27] FormSplitter: don't split coefficients --- firedrake/ffc_interface.py | 6 ------ 1 file changed, 6 deletions(-) diff --git a/firedrake/ffc_interface.py b/firedrake/ffc_interface.py index bb85a67f3f..3eee23c0c7 100644 --- a/firedrake/ffc_interface.py +++ b/firedrake/ffc_interface.py @@ -94,12 +94,6 @@ def argument(self, o): return as_vector(args) return o - def coefficient(self, o): - """Split a coefficient into its constituent spaces.""" - if isinstance(o.function_space(), types.MixedFunctionSpace): - return as_vector(list(o.split())) - return o - class FFCKernel(DiskCached): From 5dee7a4bd48eb5120cd91485b6dadf4d65745f0b Mon Sep 17 00:00:00 2001 From: Florian Rathgeber Date: Fri, 14 Feb 2014 14:22:30 +0000 Subject: [PATCH 15/27] Only sum integrands if they are on the same domain --- firedrake/ffc_interface.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/firedrake/ffc_interface.py b/firedrake/ffc_interface.py index 3eee23c0c7..4cc91701eb 100644 --- a/firedrake/ffc_interface.py +++ b/firedrake/ffc_interface.py @@ -1,6 +1,7 @@ """Provides the interface to FFC for compiling a form, and transforms the FFC- generated code in order to make it suitable for passing to the backends.""" +from collections import defaultdict from hashlib import md5 from operator import add import os @@ -47,8 +48,11 @@ def _check_version(): def sum_integrands(form): """Produce a form with the integrands on the same measure summed.""" + integrals = defaultdict(list) + for integral in form.integrals(): + integrals[integral.measure()].append(integral) return Form([it[0].reconstruct(reduce(add, [i.integrand() for i in it])) - for d, it in form.integral_groups().items()]) + for it in integrals.values()]) class FormSplitter(ReuseTransformer): From 67a39f60ba4fc95d15462ec1e7b897c23cd740fa Mon Sep 17 00:00:00 2001 From: Florian Rathgeber Date: Fri, 14 Feb 2014 14:24:44 +0000 Subject: [PATCH 16/27] Only call FormSplitter if we have a mixed form --- firedrake/ffc_interface.py | 20 +++++++++++--------- 1 file changed, 11 insertions(+), 9 deletions(-) diff --git a/firedrake/ffc_interface.py b/firedrake/ffc_interface.py index 4cc91701eb..1803b48612 100644 --- a/firedrake/ffc_interface.py +++ b/firedrake/ffc_interface.py @@ -61,12 +61,9 @@ class FormSplitter(ReuseTransformer): def split(self, form): """Split the given form.""" - fd = form.compute_form_data() - # If there is no mixed element involved, return a form per integral - if all(isinstance(e, (FiniteElement, VectorElement)) for e in fd.unique_sub_elements): - return [[((0, 0), Form([i]))] for i in sum_integrands(form).integrals()] - # Otherwise visit each integrand and obtain the tuple of sub forms - shape = tuple(len(a.function_space()) for a in fd.original_arguments) + # Visit each integrand and obtain the tuple of sub forms + shape = tuple(len(a.function_space()) + for a in form.form_data().original_arguments) forms_list = [] for it in sum_integrands(form).integrals(): forms = [] @@ -143,13 +140,18 @@ def compile_form(form, name): if not isinstance(form, Form): form = as_form(form) + fd = form.compute_form_data() + # If there is no mixed element involved, return the kernels FFC produces + if all(isinstance(e, (FiniteElement, VectorElement)) for e in fd.unique_sub_elements): + return [((0, 0), ida.integrals[0].measure(), fd.original_coefficients, kernel) + for ida, kernel in zip(fd.integral_data, FFCKernel(form, name).kernels)] + # Otherwise pre-split the form into mixed blocks before calling FFC kernels = [] for forms in FormSplitter().split(form): for (i, j), form in forms: - fd = form.compute_form_data() kernel, = FFCKernel(form, name + str(i) + str(j)).kernels - ida = fd.integral_data[0] - kernels.append(((i, j), ida.integrals[0].measure(), + fd = form.form_data() + kernels.append(((i, j), fd.integral_data[0].integrals[0].measure(), fd.original_coefficients, kernel)) return kernels From 4278503ddb1b2c3578292694dbf1445f6e3de837 Mon Sep 17 00:00:00 2001 From: Florian Rathgeber Date: Fri, 14 Mar 2014 18:27:03 +0000 Subject: [PATCH 17/27] Access integrals via preprocessed form When building integral_data for the form_data object, the domain_data which Firedrake uses to carry the coordinate field is stripped. We therefore need to get the integrals with intact integral_data via the preprocessed form. --- firedrake/ffc_interface.py | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/firedrake/ffc_interface.py b/firedrake/ffc_interface.py index 1803b48612..49a0382918 100644 --- a/firedrake/ffc_interface.py +++ b/firedrake/ffc_interface.py @@ -118,16 +118,16 @@ def __init__(self, form, name): ffc_tree = ffc_compile_form(form, prefix=name, parameters=ffc_parameters) kernels = [] - for ida, kernel in zip(form.form_data().integral_data, ffc_tree): + for it, kernel in zip(form.form_data().preprocessed_form.integrals(), ffc_tree): # Set optimization options - opts = {} if ida.domain_type not in ['cell'] else \ + opts = {} if it.domain_type() not in ['cell'] else \ {'licm': False, 'tile': None, 'vect': None, 'ap': False, 'split': None} kernels.append(Kernel(Root([incl, kernel]), '%s_%s_integral_0_%s' % - (name, ida.domain_type, ida.domain_id), opts, inc)) + (name, it.domain_type(), it.domain_id()), opts, inc)) self.kernels = tuple(kernels) self._initialized = True @@ -143,15 +143,16 @@ def compile_form(form, name): fd = form.compute_form_data() # If there is no mixed element involved, return the kernels FFC produces if all(isinstance(e, (FiniteElement, VectorElement)) for e in fd.unique_sub_elements): - return [((0, 0), ida.integrals[0].measure(), fd.original_coefficients, kernel) - for ida, kernel in zip(fd.integral_data, FFCKernel(form, name).kernels)] + return [((0, 0), it.measure(), fd.original_coefficients, kernel) + for it, kernel in zip(fd.preprocessed_form.integrals(), + FFCKernel(form, name).kernels)] # Otherwise pre-split the form into mixed blocks before calling FFC kernels = [] for forms in FormSplitter().split(form): for (i, j), form in forms: kernel, = FFCKernel(form, name + str(i) + str(j)).kernels fd = form.form_data() - kernels.append(((i, j), fd.integral_data[0].integrals[0].measure(), + kernels.append(((i, j), fd.preprocessed_form.integrals()[0].measure(), fd.original_coefficients, kernel)) return kernels From b0d25e66d179d8cdbacff470a6e4f1c2ea7abb87 Mon Sep 17 00:00:00 2001 From: Florian Rathgeber Date: Fri, 14 Mar 2014 18:27:50 +0000 Subject: [PATCH 18/27] Extract coord field from measure in assemble thunk --- firedrake/solving.py | 15 +++++---------- 1 file changed, 5 insertions(+), 10 deletions(-) diff --git a/firedrake/solving.py b/firedrake/solving.py index 6e9b228b16..3d235cb60a 100644 --- a/firedrake/solving.py +++ b/firedrake/solving.py @@ -403,9 +403,6 @@ def _assemble(f, tensor=None, bcs=None): integrals = fd.preprocessed_form.integrals() - # Extract coordinate field - coords = integrals[0].measure().domain_data() - def get_rank(arg): return arg.function_space().rank has_vec_fs = lambda arg: isinstance(arg.function_space(), types.VectorFunctionSpace) @@ -423,7 +420,6 @@ def mixed_plus_vfs_error(arg): mixed_plus_vfs_error(test) mixed_plus_vfs_error(trial) - m = test.function_space().mesh() map_pairs = [] cell_domains = [] exterior_facet_domains = [] @@ -491,7 +487,6 @@ def mixed_plus_vfs_error(arg): elif is_vec: test = fd.original_arguments[0] mixed_plus_vfs_error(test) - m = test.function_space().mesh() if tensor is None: result_function = types.Function(test.function_space()) tensor = result_function.dat @@ -501,8 +496,6 @@ def mixed_plus_vfs_error(arg): tensor.zero() result = lambda: result_function else: - m = coords.function_space().mesh() - # 0-forms are always scalar if tensor is None: tensor = op2.Global(1, [0.0]) @@ -526,6 +519,8 @@ def thunk(bcs): top = any(bc.sub_domain == "top" for bc in bcs) extruded_bcs = (bottom, top) for (i, j), measure, coefficients, kernel in kernels: + coords = measure.domain_data() + m = coords.function_space().mesh() # Extract block from tensor and test/trial spaces # FIXME Ugly variable renaming required because functions are not # lexical closures in Python and we're writing to these variables @@ -623,7 +618,7 @@ def thunk(bcs): #domain id: interior horizontal, bottom or top. #Get the list of sets and globals required for parallel loop construction. - set_global_list = m.exterior_facets.measure_set(integral.measure()) + set_global_list = m.exterior_facets.measure_set(measure) #Iterate over the list and assemble all the args of the parallel loop for (index, set) in set_global_list: @@ -656,7 +651,7 @@ def thunk(bcs): else: tensor_arg = t(op2.INC) - args = [kernel, m.exterior_facets.measure_set(integral.measure()), tensor_arg, + args = [kernel, m.exterior_facets.measure_set(measure), tensor_arg, coords.dat(op2.READ, coords.exterior_facet_node_map(), flatten=True)] for c in coefficients: @@ -712,7 +707,7 @@ def thunk(bcs): else: tensor_arg = t(op2.INC) - args = [kernel, m.interior_facets.measure_set(integral.measure()), tensor_arg, + args = [kernel, m.interior_facets.measure_set(measure), tensor_arg, coords.dat(op2.READ, coords.cell_node_map(), flatten=True)] for c in coefficients: From 6ad536f5e385cc0f8a871190439a7b34ff3c4a55 Mon Sep 17 00:00:00 2001 From: Florian Rathgeber Date: Fri, 14 Mar 2014 18:32:01 +0000 Subject: [PATCH 19/27] Add firedrake-clean script to remove cached FFC kernels --- scripts/firedrake-clean | 8 ++++++++ {tools => scripts}/gmsh2triangle.py | 0 setup.py | 2 ++ 3 files changed, 10 insertions(+) create mode 100755 scripts/firedrake-clean rename {tools => scripts}/gmsh2triangle.py (100%) diff --git a/scripts/firedrake-clean b/scripts/firedrake-clean new file mode 100755 index 0000000000..785d55cc9c --- /dev/null +++ b/scripts/firedrake-clean @@ -0,0 +1,8 @@ +#!/usr/bin/env python + +from firedrake.ffc_interface import clear_cache, FFCKernel + + +if __name__ == '__main__': + print 'Removing cached ffc kernels from %s' % FFCKernel._cachedir + clear_cache() diff --git a/tools/gmsh2triangle.py b/scripts/gmsh2triangle.py similarity index 100% rename from tools/gmsh2triangle.py rename to scripts/gmsh2triangle.py diff --git a/setup.py b/setup.py index a175db9a04..8b9c3a3dab 100644 --- a/setup.py +++ b/setup.py @@ -1,5 +1,6 @@ from distutils.core import setup from distutils.extension import Extension +from glob import glob import numpy as np try: @@ -25,6 +26,7 @@ author_email="firedrake@imperial.ac.uk", url="http://firedrakeproject.org", packages=["firedrake", "evtk"], + scripts=glob('scripts/*'), ext_modules=[Extension('firedrake.core_types', sources=firedrake_sources, include_dirs=[np.get_include()]), From 6de929b5296b374335431719c6d9f400acdc0d35 Mon Sep 17 00:00:00 2001 From: Florian Rathgeber Date: Mon, 17 Mar 2014 17:45:32 +0000 Subject: [PATCH 20/27] Form splitter: correctly deal with renumbered form arguments --- firedrake/ffc_interface.py | 30 ++++++++++++++++++++++-------- 1 file changed, 22 insertions(+), 8 deletions(-) diff --git a/firedrake/ffc_interface.py b/firedrake/ffc_interface.py index 49a0382918..ab4a3aa218 100644 --- a/firedrake/ffc_interface.py +++ b/firedrake/ffc_interface.py @@ -62,17 +62,31 @@ class FormSplitter(ReuseTransformer): def split(self, form): """Split the given form.""" # Visit each integrand and obtain the tuple of sub forms - shape = tuple(len(a.function_space()) - for a in form.form_data().original_arguments) + args = tuple((a.count(), len(a.function_space())) + for a in form.form_data().original_arguments) forms_list = [] for it in sum_integrands(form).integrals(): forms = [] - for i in range(shape[0] if len(shape) > 0 else 1): - for j in range(shape[1] if len(shape) > 1 else 1): - self._idx = {-2: i, -1: j} - integrand = self.visit(it.integrand()) - if not isinstance(integrand, Zero): - forms.append([((i, j), integrand * it.measure())]) + + def visit(idx): + integrand = self.visit(it.integrand()) + if not isinstance(integrand, Zero): + forms.append([(idx, integrand * it.measure())]) + # 0 form + if not args: + visit((0, 0)) + # 1 form + elif len(args) == 1: + count, l = args[0] + for i in range(l): + self._idx = {count: i} + visit((i, 0)) + # 2 form + elif len(args) == 2: + for i in range(args[0][1]): + for j in range(args[1][1]): + self._idx = {args[0][0]: i, args[1][0]: j} + visit((i, j)) forms_list += forms return forms_list From 92d066fb2efc6e434bb11b4fda6f61b35445c278 Mon Sep 17 00:00:00 2001 From: Florian Rathgeber Date: Mon, 17 Mar 2014 17:52:34 +0000 Subject: [PATCH 21/27] Add mixed helmholtz test cast in nonlinear form --- tests/regression/test_helmholtz_mixed.py | 20 +++++++++++++------- 1 file changed, 13 insertions(+), 7 deletions(-) diff --git a/tests/regression/test_helmholtz_mixed.py b/tests/regression/test_helmholtz_mixed.py index acb60c7fec..5ae1a732b1 100644 --- a/tests/regression/test_helmholtz_mixed.py +++ b/tests/regression/test_helmholtz_mixed.py @@ -2,7 +2,7 @@ from firedrake import * -def helmholtz_mixed(x, V1, V2): +def helmholtz_mixed(x, V1, V2, action=False): # Create mesh and define function space mesh = UnitSquareMesh(2**x, 2**x) V1 = FunctionSpace(mesh, *V1, name="V") @@ -22,11 +22,16 @@ def helmholtz_mixed(x, V1, V2): # Compute solution x = Function(W) + if action: + system = action(a, x) - L == 0 + else: + system = a == L + # Block system is: # V Ct # Ch P # Eliminate V by forming a schur complement - solve(a == L, x, solver_parameters={'pc_type': 'fieldsplit', + solve(system, x, solver_parameters={'pc_type': 'fieldsplit', 'pc_fieldsplit_type': 'schur', 'ksp_type': 'cg', 'pc_fieldsplit_schur_fact_type': 'FULL', @@ -38,11 +43,12 @@ def helmholtz_mixed(x, V1, V2): return sqrt(assemble(dot(x[2] - f, x[2] - f) * dx)) -@pytest.mark.parametrize(('V1', 'V2', 'threshold'), - [(('RT', 1), ('DG', 0), 1.9), - (('BDM', 1), ('DG', 0), 1.89), - (('BDFM', 2), ('DG', 1), 1.9)]) -def test_firedrake_helmholtz(V1, V2, threshold): +@pytest.mark.parametrize(('V1', 'V2', 'threshold', 'action'), + [(('RT', 1), ('DG', 0), 1.9, False), + (('BDM', 1), ('DG', 0), 1.89, False), + (('BDM', 1), ('DG', 0), 1.89, True), + (('BDFM', 2), ('DG', 1), 1.9, False)]) +def test_firedrake_helmholtz(V1, V2, threshold, action): import numpy as np diff = np.array([helmholtz_mixed(i, V1, V2) for i in range(3, 6)]) print "l2 error norms:", diff From 79561227cbc82700f2733ee7dfad0044e90cbcb9 Mon Sep 17 00:00:00 2001 From: Florian Rathgeber Date: Fri, 21 Mar 2014 15:08:03 +0000 Subject: [PATCH 22/27] Add test assembling derivative of form with 0 block --- tests/regression/test_split.py | 22 ++++++++++++++++++++++ 1 file changed, 22 insertions(+) create mode 100644 tests/regression/test_split.py diff --git a/tests/regression/test_split.py b/tests/regression/test_split.py new file mode 100644 index 0000000000..209b5720dc --- /dev/null +++ b/tests/regression/test_split.py @@ -0,0 +1,22 @@ +import pytest +from firedrake import * + + +def test_assemble_split_derivative(): + """Assemble the derivative of a form with a zero block.""" + mesh = UnitSquareMesh(1, 1) + V1 = FunctionSpace(mesh, "BDM", 1, name="V") + V2 = FunctionSpace(mesh, "DG", 0, name="P") + W = V1 * V2 + + x = Function(W) + u, p = split(x) + v, q = TestFunctions(W) + + F = (inner(u, v) + v[1]*p)*dx + + assert assemble(derivative(F, x)) + +if __name__ == '__main__': + import os + pytest.main(os.path.abspath(__file__)) From 92728319a5415705d396413817ddeb8fd40c9420 Mon Sep 17 00:00:00 2001 From: Florian Rathgeber Date: Tue, 1 Apr 2014 15:41:33 +0100 Subject: [PATCH 23/27] Cache split arguments to keep them unique --- firedrake/ffc_interface.py | 20 +++++++++++++------- 1 file changed, 13 insertions(+), 7 deletions(-) diff --git a/firedrake/ffc_interface.py b/firedrake/ffc_interface.py index ab4a3aa218..d588bb6c8f 100644 --- a/firedrake/ffc_interface.py +++ b/firedrake/ffc_interface.py @@ -80,34 +80,40 @@ def visit(idx): count, l = args[0] for i in range(l): self._idx = {count: i} + self._args = {} visit((i, 0)) # 2 form elif len(args) == 2: for i in range(args[0][1]): for j in range(args[1][1]): self._idx = {args[0][0]: i, args[1][0]: j} + self._args = {} visit((i, j)) forms_list += forms return forms_list - def argument(self, o): + def argument(self, arg): """Split an argument into its constituent spaces.""" - if isinstance(o.function_space(), types.MixedFunctionSpace): + if isinstance(arg.function_space(), types.MixedFunctionSpace): + if arg in self._args: + return self._args[arg] args = [] - for i, fs in enumerate(o.function_space().split()): - a = Argument(fs.ufl_element(), fs, o.count()) + for i, fs in enumerate(arg.function_space().split()): + # Look up the split argument in cache since we want it unique + a = Argument(fs.ufl_element(), fs, arg.count()) if a.shape(): - if self._idx[o.count()] == i: + if self._idx[arg.count()] == i: args += [a[j] for j in range(a.shape()[0])] else: args += [Zero() for j in range(a.shape()[0])] else: - if self._idx[o.count()] == i: + if self._idx[arg.count()] == i: args.append(a) else: args.append(Zero()) + self._args[arg] = as_vector(args) return as_vector(args) - return o + return arg class FFCKernel(DiskCached): From e16e4dd0448694ac9ea0090df73f4b4c7773eebd Mon Sep 17 00:00:00 2001 From: Florian Rathgeber Date: Wed, 2 Apr 2014 16:38:29 +0100 Subject: [PATCH 24/27] Allow extracting subspace 0 from (Vector)FunctionSpace --- firedrake/types.py | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/firedrake/types.py b/firedrake/types.py index e1cfd6c1f9..a6d424bc4c 100644 --- a/firedrake/types.py +++ b/firedrake/types.py @@ -98,6 +98,11 @@ def _process_args(cls, *args, **kwargs): def _cache_key(cls, mesh, family, degree=None, name=None, vfamily=None, vdegree=None): return family, degree, vfamily, vdegree + def __getitem__(self, i): + """Return self if ``i`` is 0, otherwise raise an error.""" + assert i == 0, "Can only extract subspace 0 from %r" % self + return self + class VectorFunctionSpace(FunctionSpaceBase): """A vector finite element :class:`FunctionSpace`.""" @@ -135,6 +140,11 @@ def _process_args(cls, *args, **kwargs): def _cache_key(cls, mesh, family, degree=None, dim=None, name=None, vfamily=None, vdegree=None): return family, degree, dim, vfamily, vdegree + def __getitem__(self, i): + """Return self if ``i`` is 0, otherwise raise an error.""" + assert i == 0, "Can only extract subspace 0 from %r" % self + return self + class MixedFunctionSpace(FunctionSpaceBase): """A mixed finite element :class:`FunctionSpace`.""" From e1ffde0f7f27eb8b1911f021998dca3589a334d4 Mon Sep 17 00:00:00 2001 From: Florian Rathgeber Date: Wed, 2 Apr 2014 16:45:52 +0100 Subject: [PATCH 25/27] Add helper functions to build mat/vec in assemble thunk --- firedrake/solving.py | 107 +++++++++++++++++++------------------------ 1 file changed, 47 insertions(+), 60 deletions(-) diff --git a/firedrake/solving.py b/firedrake/solving.py index 3d235cb60a..cf5131aa20 100644 --- a/firedrake/solving.py +++ b/firedrake/solving.py @@ -483,6 +483,12 @@ def mixed_plus_vfs_error(arg): result_matrix.bcs = bcs tensor = tensor._M tensor.zero() + + def mat(testmap, trialmap, i, j): + return tensor[i, j](op2.INC, + (testmap(test.function_space()[i])[op2.i[0]], + trialmap(trial.function_space()[j])[op2.i[1]]), + flatten=has_vec_fs(test)) result = lambda: result_matrix elif is_vec: test = fd.original_arguments[0] @@ -493,7 +499,12 @@ def mixed_plus_vfs_error(arg): else: result_function = tensor tensor = result_function.dat - tensor.zero() + tensor.zero() + + def vec(testmap, i): + return tensor[i](op2.INC, + testmap(test.function_space()[i])[op2.i[0]], + flatten=has_vec_fs(test)) result = lambda: result_function else: # 0-forms are always scalar @@ -525,33 +536,19 @@ def thunk(bcs): # FIXME Ugly variable renaming required because functions are not # lexical closures in Python and we're writing to these variables if is_mat and tensor.sparsity.shape > (1, 1): - t = tensor[i, j] - ts = test.function_space()[i] - tr = trial.function_space()[j] tsbc = [bc for bc in bcs if bc.function_space().index == i] trbc = [bc for bc in bcs if bc.function_space().index == j] elif is_mat: - t = tensor - ts, tr = test, trial tsbc, trbc = bcs, bcs - elif is_vec and len(tensor) > 1: - t = tensor[i] - ts = test.function_space()[i] - elif is_vec: - t = tensor - ts = test - else: - t = tensor if measure.domain_type() == 'cell': if is_mat: - tensor_arg = t(op2.INC, (ts.cell_node_map(tsbc)[op2.i[0]], - tr.cell_node_map(trbc)[op2.i[1]]), - flatten=has_vec_fs(test)) + tensor_arg = mat(lambda s: s.cell_node_map(tsbc), + lambda s: s.cell_node_map(trbc), + i, j) elif is_vec: - tensor_arg = t(op2.INC, ts.cell_node_map()[op2.i[0]], - flatten=has_vec_fs(test)) + tensor_arg = vec(lambda s: s.cell_node_map(), i) else: - tensor_arg = t(op2.INC) + tensor_arg = tensor(op2.INC) itspace = m.cell_set itspace._extruded_bcs = extruded_bcs @@ -569,6 +566,7 @@ def thunk(bcs): op2.par_loop(*args) except MapValueError: raise RuntimeError("Integral measure does not match measure of all coefficients/arguments") + elif measure.domain_type() == 'exterior_facet': if op2.MPI.parallel: raise \ @@ -576,15 +574,13 @@ def thunk(bcs): "No support for facet integrals under MPI yet") if is_mat: - tensor_arg = t(op2.INC, (ts.exterior_facet_node_map(tsbc)[op2.i[0]], - tr.exterior_facet_node_map(trbc)[op2.i[1]]), - flatten=has_vec_fs(test)) + tensor_arg = mat(lambda s: s.exterior_facet_node_map(tsbc), + lambda s: s.exterior_facet_node_map(trbc), + i, j) elif is_vec: - tensor_arg = t(op2.INC, - ts.exterior_facet_node_map()[op2.i[0]], - flatten=has_vec_fs(test)) + tensor_arg = vec(lambda s: s.exterior_facet_node_map(), i) else: - tensor_arg = t(op2.INC) + tensor_arg = tensor(op2.INC) args = [kernel, m.exterior_facets.measure_set(measure), tensor_arg, coords.dat(op2.READ, coords.exterior_facet_node_map(), flatten=True)] @@ -604,12 +600,11 @@ def thunk(bcs): "No support for facet integrals under MPI yet") if is_mat: - tensor_arg = t(op2.INC, (ts.cell_node_map(tsbc)[op2.i[0]], - tr.cell_node_map(trbc)[op2.i[1]]), - flatten=has_vec_fs(test)) + tensor_arg = mat(lambda s: s.cell_node_map(tsbc), + lambda s: s.cell_node_map(trbc), + i, j) elif is_vec: - tensor_arg = t(op2.INC, ts.cell_node_map()[op2.i[0]], - flatten=has_vec_fs(test)) + tensor_arg = vec(lambda s: s.cell_node_map(), i) else: tensor_arg = t(op2.INC) @@ -641,13 +636,11 @@ def thunk(bcs): "No support for facet integrals under MPI yet") if is_mat: - tensor_arg = t(op2.INC, - (ts.exterior_facet_node_map(tsbc)[op2.i[0]], - tr.exterior_facet_node_map(trbc)[op2.i[1]]), - flatten=has_vec_fs(test)) + tensor_arg = mat(lambda s: s.exterior_facet_node_map(tsbc), + lambda s: s.exterior_facet_node_map(trbc), + i, j) elif is_vec: - tensor_arg = t(op2.INC, ts.exterior_facet_node_map()[op2.i[0]], - flatten=has_vec_fs(test)) + tensor_arg = vec(lambda s: s.exterior_facet_node_map(), i) else: tensor_arg = t(op2.INC) @@ -670,15 +663,13 @@ def thunk(bcs): "No support for facet integrals under MPI yet") if is_mat: - tensor_arg = t(op2.INC, - (ts.interior_facet_node_map(tsbc)[op2.i[0]], - tr.interior_facet_node_map(trbc)[op2.i[1]]), - flatten=True) + tensor_arg = mat(lambda s: s.interior_facet_node_map(tsbc), + lambda s: s.interior_facet_node_map(trbc), + i, j) elif is_vec: - tensor_arg = t(op2.INC, ts.interior_facet_node_map()[op2.i[0]], - flatten=True) + tensor_arg = vec(lambda s: s.interior_facet_node_map(), i) else: - tensor_arg = t(op2.INC) + tensor_arg = tensor(op2.INC) args = [kernel, m.interior_facets.set, tensor_arg, coords.dat(op2.READ, coords.interior_facet_node_map(), flatten=True)] @@ -698,14 +689,13 @@ def thunk(bcs): "No support for facet integrals under MPI yet") if is_mat: - tensor_arg = t(op2.INC, (ts.cell_node_map(tsbc)[op2.i[0]], - tr.cell_node_map(trbc)[op2.i[1]]), - flatten=True) + tensor_arg = mat(lambda s: s.cell_node_map(tsbc), + lambda s: s.cell_node_map(trbc), + i, j) elif is_vec: - tensor_arg = t(op2.INC, ts.cell_node_map()[op2.i[0]], - flatten=True) + tensor_arg = vec(lambda s: s.cell_node_map(), i) else: - tensor_arg = t(op2.INC) + tensor_arg = tensor(op2.INC) args = [kernel, m.interior_facets.measure_set(measure), tensor_arg, coords.dat(op2.READ, coords.cell_node_map(), @@ -725,16 +715,13 @@ def thunk(bcs): "No support for facet integrals under MPI yet") if is_mat: - tensor_arg = t(op2.INC, - (ts.interior_facet_node_map(tsbc)[op2.i[0]], - tr.interior_facet_node_map(trbc)[op2.i[1]]), - flatten=True) + tensor_arg = mat(lambda s: s.interior_facet_node_map(tsbc), + lambda s: s.interior_facet_node_map(trbc), + i, j) elif is_vec: - tensor_arg = t(op2.INC, - ts.interior_facet_node_map()[op2.i[0]], - flatten=True) + tensor_arg = vec(lambda s: s.interior_facet_node_map(), i) else: - tensor_arg = t(op2.INC) + tensor_arg = tensor(op2.INC) args = [kernel, m.interior_facets.set, tensor_arg, coords.dat(op2.READ, coords.interior_facet_node_map(), flatten=True)] @@ -759,7 +746,7 @@ def thunk(bcs): # block is on the matrix diagonal and its index matches the # index of the function space the bc is defined on. if i == j and (fs.index is None or fs.index == i): - t.inc_local_diagonal_entries(bc.nodes) + tensor[i, j].inc_local_diagonal_entries(bc.nodes) return result() From b202a7b79adcfbb39beb1e5594c771aeb695075b Mon Sep 17 00:00:00 2001 From: Florian Rathgeber Date: Wed, 2 Apr 2014 18:04:35 +0100 Subject: [PATCH 26/27] Eliminate duplicate code for horizontal interior/exterior facet --- firedrake/solving.py | 59 +++----------------------------------------- 1 file changed, 3 insertions(+), 56 deletions(-) diff --git a/firedrake/solving.py b/firedrake/solving.py index cf5131aa20..8bb429ce7f 100644 --- a/firedrake/solving.py +++ b/firedrake/solving.py @@ -567,7 +567,7 @@ def thunk(bcs): except MapValueError: raise RuntimeError("Integral measure does not match measure of all coefficients/arguments") - elif measure.domain_type() == 'exterior_facet': + elif measure.domain_type() in ['exterior_facet', 'exterior_facet_vert']: if op2.MPI.parallel: raise \ NotImplementedError( @@ -606,7 +606,7 @@ def thunk(bcs): elif is_vec: tensor_arg = vec(lambda s: s.cell_node_map(), i) else: - tensor_arg = t(op2.INC) + tensor_arg = tensor(op2.INC) #In the case of extruded meshes with horizontal facet integrals, two #parallel loops will (potentially) get created and called based on the @@ -629,34 +629,7 @@ def thunk(bcs): except MapValueError: raise RuntimeError("Integral measure does not match measure of all coefficients/arguments") - elif measure.domain_type() == 'exterior_facet_vert': - if op2.MPI.parallel: - raise \ - NotImplementedError( - "No support for facet integrals under MPI yet") - - if is_mat: - tensor_arg = mat(lambda s: s.exterior_facet_node_map(tsbc), - lambda s: s.exterior_facet_node_map(trbc), - i, j) - elif is_vec: - tensor_arg = vec(lambda s: s.exterior_facet_node_map(), i) - else: - tensor_arg = t(op2.INC) - - args = [kernel, m.exterior_facets.measure_set(measure), tensor_arg, - coords.dat(op2.READ, coords.exterior_facet_node_map(), - flatten=True)] - for c in coefficients: - args.append(c.dat(op2.READ, c.exterior_facet_node_map(), - flatten=has_vec_fs(c))) - args.append(m.exterior_facets.local_facet_dat(op2.READ)) - try: - op2.par_loop(*args) - except MapValueError: - raise RuntimeError("Integral measure does not match measure of all coefficients/arguments") - - elif measure.domain_type() == 'interior_facet': + elif measure.domain_type() in ['interior_facet', 'interior_facet_vert']: if op2.MPI.parallel: raise \ NotImplementedError( @@ -708,32 +681,6 @@ def thunk(bcs): except MapValueError: raise RuntimeError("Integral measure does not match measure of all coefficients/arguments") - elif measure.domain_type() == 'interior_facet_vert': - if op2.MPI.parallel: - raise \ - NotImplementedError( - "No support for facet integrals under MPI yet") - - if is_mat: - tensor_arg = mat(lambda s: s.interior_facet_node_map(tsbc), - lambda s: s.interior_facet_node_map(trbc), - i, j) - elif is_vec: - tensor_arg = vec(lambda s: s.interior_facet_node_map(), i) - else: - tensor_arg = tensor(op2.INC) - args = [kernel, m.interior_facets.set, tensor_arg, - coords.dat(op2.READ, coords.interior_facet_node_map(), - flatten=True)] - for c in coefficients: - args.append(c.dat(op2.READ, c.interior_facet_node_map(), - flatten=has_vec_fs(c))) - args.append(m.interior_facets.local_facet_dat(op2.READ)) - try: - op2.par_loop(*args) - except MapValueError: - raise RuntimeError("Integral measure does not match measure of all coefficients/arguments") - else: raise RuntimeError('Unknown domain type "%s"' % measure.domain_type()) From a15870451ae0e7dcf9cab014418e1bcbfcf6cbd4 Mon Sep 17 00:00:00 2001 From: Florian Rathgeber Date: Wed, 2 Apr 2014 18:28:03 +0100 Subject: [PATCH 27/27] Allow overriding FFC kernel cache directory from environment --- firedrake/ffc_interface.py | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/firedrake/ffc_interface.py b/firedrake/ffc_interface.py index d588bb6c8f..ce42c29733 100644 --- a/firedrake/ffc_interface.py +++ b/firedrake/ffc_interface.py @@ -4,7 +4,7 @@ from collections import defaultdict from hashlib import md5 from operator import add -import os +from os import path, environ, getuid, makedirs import tempfile from ufl import Form, FiniteElement, VectorElement, as_vector @@ -30,7 +30,7 @@ ffc_parameters['pyop2-ir'] = True # Include an md5 hash of firedrake_geometry.h in the cache key -with open(os.path.join(os.path.dirname(__file__), 'firedrake_geometry.h')) as f: +with open(path.join(path.dirname(__file__), 'firedrake_geometry.h')) as f: _firedrake_geometry_md5 = md5(f.read()).hexdigest() @@ -119,8 +119,9 @@ def argument(self, arg): class FFCKernel(DiskCached): _cache = {} - _cachedir = os.path.join(tempfile.gettempdir(), - 'firedrake-ffc-kernel-cache-uid%d' % os.getuid()) + _cachedir = environ.get('FIREDRAKE_FFC_KERNEL_CACHE_DIR', + path.join(tempfile.gettempdir(), + 'firedrake-ffc-kernel-cache-uid%d' % getuid())) @classmethod def _cache_key(cls, form, name): @@ -134,7 +135,7 @@ def __init__(self, form, name): return incl = PreprocessNode('#include "firedrake_geometry.h"\n') - inc = [os.path.dirname(__file__)] + inc = [path.dirname(__file__)] ffc_tree = ffc_compile_form(form, prefix=name, parameters=ffc_parameters) kernels = [] @@ -181,7 +182,7 @@ def clear_cache(): """Clear the PyOP2 FFC kernel cache.""" if MPI.comm.rank != 0: return - if os.path.exists(FFCKernel._cachedir): + if path.exists(FFCKernel._cachedir): import shutil shutil.rmtree(FFCKernel._cachedir, ignore_errors=True) _ensure_cachedir() @@ -189,8 +190,8 @@ def clear_cache(): def _ensure_cachedir(): """Ensure that the FFC kernel cache directory exists.""" - if not os.path.exists(FFCKernel._cachedir) and MPI.comm.rank == 0: - os.makedirs(FFCKernel._cachedir) + if not path.exists(FFCKernel._cachedir) and MPI.comm.rank == 0: + makedirs(FFCKernel._cachedir) _check_version() _ensure_cachedir()