diff --git a/firedrake/ffc_interface.py b/firedrake/ffc_interface.py new file mode 100644 index 0000000000..ce42c29733 --- /dev/null +++ b/firedrake/ffc_interface.py @@ -0,0 +1,197 @@ +"""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 +from os import path, environ, getuid, makedirs +import tempfile + +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 +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 + +import types + +_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(path.join(path.dirname(__file__), 'firedrake_geometry.h')) as f: + _firedrake_geometry_md5 = md5(f.read()).hexdigest() + + +def _check_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]: + return + except AttributeError: + pass + raise RuntimeError("Incompatible PyOP2 version %s and FFC PyOP2 version %s." + % (version, getattr(constants, 'PYOP2_VERSION', 'unknown'))) + + +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 it in integrals.values()]) + + +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.""" + # Visit each integrand and obtain the tuple of sub forms + 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 = [] + + 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} + 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, arg): + """Split an argument into its constituent spaces.""" + if isinstance(arg.function_space(), types.MixedFunctionSpace): + if arg in self._args: + return self._args[arg] + args = [] + 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[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[arg.count()] == i: + args.append(a) + else: + args.append(Zero()) + self._args[arg] = as_vector(args) + return as_vector(args) + return arg + + +class FFCKernel(DiskCached): + + _cache = {} + _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): + 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 = [path.dirname(__file__)] + ffc_tree = ffc_compile_form(form, prefix=name, parameters=ffc_parameters) + + kernels = [] + for it, kernel in zip(form.form_data().preprocessed_form.integrals(), ffc_tree): + # Set optimization options + 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, it.domain_type(), it.domain_id()), opts, inc)) + self.kernels = tuple(kernels) + self._initialized = True + + +def compile_form(form, name): + """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): + 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), 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.preprocessed_form.integrals()[0].measure(), + fd.original_coefficients, kernel)) + return kernels + + +def clear_cache(): + """Clear the PyOP2 FFC kernel cache.""" + if MPI.comm.rank != 0: + return + if 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 path.exists(FFCKernel._cachedir) and MPI.comm.rank == 0: + 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..8bb429ce7f 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() @@ -402,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) @@ -422,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 = [] @@ -486,22 +483,30 @@ 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] 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 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: - m = coords.function_space().mesh() - # 0-forms are always scalar if tensor is None: tensor = op2.Global(1, [0.0]) @@ -524,16 +529,24 @@ 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: + 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 + if is_mat and tensor.sparsity.shape > (1, 1): + 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: + tsbc, trbc = bcs, bcs + 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 = mat(lambda s: s.cell_node_map(tsbc), + lambda s: s.cell_node_map(trbc), + i, j) elif is_vec: - tensor_arg = tensor(op2.INC, test.cell_node_map()[op2.i[0]], - flatten=has_vec_fs(test)) + tensor_arg = vec(lambda s: s.cell_node_map(), i) else: tensor_arg = tensor(op2.INC) @@ -545,7 +558,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))) @@ -553,27 +566,25 @@ 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() in ['exterior_facet', '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 = mat(lambda s: s.exterior_facet_node_map(tsbc), + lambda s: s.exterior_facet_node_map(trbc), + i, j) elif is_vec: - tensor_arg = tensor(op2.INC, - test.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 = tensor(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 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)) @@ -582,21 +593,18 @@ 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 = mat(lambda s: s.cell_node_map(tsbc), + lambda s: s.cell_node_map(trbc), + i, j) elif is_vec: - tensor_arg = tensor(op2.INC, - test.cell_node_map()[op2.i[0]], - flatten=has_vec_fs(test)) + tensor_arg = vec(lambda s: s.cell_node_map(), i) else: tensor_arg = tensor(op2.INC) @@ -605,7 +613,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: @@ -613,7 +621,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: @@ -621,58 +629,24 @@ def thunk(bcs): except MapValueError: raise RuntimeError("Integral measure does not match measure of all coefficients/arguments") - elif 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)) - elif is_vec: - tensor_arg = tensor(op2.INC, - test.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, - coords.dat(op2.READ, coords.exterior_facet_node_map(), - flatten=True)] - for c in fd.original_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 domain_type == 'interior_facet': + elif measure.domain_type() in ['interior_facet', '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 = mat(lambda s: s.interior_facet_node_map(tsbc), + lambda s: s.interior_facet_node_map(trbc), + i, j) elif is_vec: - tensor_arg = tensor( - op2.INC, test.interior_facet_node_map()[op2.i[0]], - flatten=True) + 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 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)) @@ -681,28 +655,25 @@ 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 = mat(lambda s: s.cell_node_map(tsbc), + lambda s: s.cell_node_map(trbc), + i, j) elif is_vec: - tensor_arg = tensor( - op2.INC, test.cell_node_map()[op2.i[0]], - flatten=True) + tensor_arg = vec(lambda s: s.cell_node_map(), i) else: tensor_arg = tensor(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 fd.original_coefficients: + for c in coefficients: args.append(c.dat(op2.READ, c.cell_node_map(), flatten=has_vec_fs(c))) try: @@ -710,51 +681,19 @@ def thunk(bcs): except MapValueError: raise RuntimeError("Integral measure does not match measure of all coefficients/arguments") - elif 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) - elif is_vec: - tensor_arg = tensor( - op2.INC, test.interior_facet_node_map()[op2.i[0]], - flatten=True) - 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 fd.original_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"' % 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): + tensor[i, j].inc_local_diagonal_entries(bc.nodes) return result() 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`.""" 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(): 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()]), 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 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__)) diff --git a/tests/test_ffc_interface.py b/tests/test_ffc_interface.py new file mode 100644 index 0000000000..60b46f66d6 --- /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[-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') + + 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') + + 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][-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][-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][-1].code and \ + 'exterior_facet_integral' in k[1][-1].code and len(k) == 2 + +if __name__ == '__main__': + pytest.main(os.path.abspath(__file__))