From c634cf96c13185486c5510677831eb4e614af297 Mon Sep 17 00:00:00 2001 From: Florian Rathgeber Date: Fri, 31 Jan 2014 17:39:34 +0000 Subject: [PATCH 01/18] Remove burgers demo This demo is broken and not used as a regression test. Prevent further bitrot by removing it. --- demo/burgers.py | 200 ------------------------------------------------ 1 file changed, 200 deletions(-) delete mode 100644 demo/burgers.py diff --git a/demo/burgers.py b/demo/burgers.py deleted file mode 100644 index c9754540f..000000000 --- a/demo/burgers.py +++ /dev/null @@ -1,200 +0,0 @@ -# This file is part of PyOP2 -# -# PyOP2 is Copyright (c) 2012, Imperial College London and -# others. Please see the AUTHORS file in the main source directory for -# a full list of copyright holders. All rights reserved. -# -# Redistribution and use in source and binary forms, with or without -# modification, are permitted provided that the following conditions -# are met: -# -# * Redistributions of source code must retain the above copyright -# notice, this list of conditions and the following disclaimer. -# * Redistributions in binary form must reproduce the above copyright -# notice, this list of conditions and the following disclaimer in the -# documentation and/or other materials provided with the distribution. -# * The name of Imperial College London or that of other -# contributors may not be used to endorse or promote products -# derived from this software without specific prior written -# permission. -# -# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTERS -# ''AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT -# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS -# FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE -# COPYRIGHT HOLDERS OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, -# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES -# (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR -# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) -# HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, -# STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) -# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED -# OF THE POSSIBILITY OF SUCH DAMAGE. - -"""Burgers equation demo (unstable forward-Euler integration) - -This demo solves the steady-state Burgers equation on a unit interval. -""" - -from pyop2 import op2, utils -from pyop2.ffc_interface import compile_form -from ufl import * -import numpy as np -import pylab - -parser = utils.parser(group=True, - description=__doc__) -parser.add_argument('-p', '--plot', - action='store_true', - help='Plot the resulting L2 error norm') - -opt = vars(parser.parse_args()) -op2.init(**opt) - -# Simulation parameters -n = 100 -nu = 0.0001 -timestep = 1.0 / n - -# Create simulation data structures - -nodes = op2.Set(n, "nodes") -b_nodes = op2.Set(2, "b_nodes") -elements = op2.Set(n - 1, "elements") - -elem_node_map = [item for sublist in [(x, x + 1) - for x in xrange(n - 1)] for item in sublist] - -elem_node = op2.Map(elements, nodes, 2, elem_node_map, "elem_node") - -b_node_node_map = [0, n - 1] -b_node_node = op2.Map(b_nodes, nodes, 1, b_node_node_map, "b_node_node") - -coord_vals = [i * (1.0 / (n - 1)) for i in xrange(n)] -coords = op2.Dat(nodes, coord_vals, np.float64, "coords") - -tracer_vals = np.asarray([0.0] * n, dtype=np.float64) -tracer = op2.Dat(nodes, tracer_vals, np.float64, "tracer") - -tracer_old_vals = np.asarray([0.0] * n, dtype=np.float64) -tracer_old = op2.Dat(nodes, tracer_old_vals, np.float64, "tracer_old") - -b_vals = np.asarray([0.0] * n, dtype=np.float64) -b = op2.Dat(nodes, b_vals, np.float64, "b") - -bdry_vals = [0.0, 1.0] -bdry = op2.Dat(nodes, bdry_vals, np.float64, "bdry") - -sparsity = op2.Sparsity((nodes, nodes), (elem_node, elem_node), "sparsity") -mat = op2.Mat(sparsity, np.float64, "mat") - -# Set up finite element problem - -V = FiniteElement("Lagrange", "interval", 1) -u = Coefficient(V) -u_next = TrialFunction(V) -v = TestFunction(V) - -a = (dot(u, grad(u_next)) * v + nu * grad(u_next) * grad(v)) * dx -L = v * u * dx - -burgers, = compile_form(a, "burgers") -rhs, = compile_form(L, "rhs") - -# Initial condition - -i_cond_code = """ -void i_cond(double *c, double *t) -{ - double pi = 3.14159265358979; - *t = *c*2; -} -""" - -i_cond = op2.Kernel(i_cond_code, "i_cond") - -op2.par_loop(i_cond, nodes, - coords(op2.READ), - tracer(op2.WRITE)) - -# Boundary condition - -strongbc_rhs = op2.Kernel( - "void strongbc_rhs(double *v, double *t) { *t = *v; }", "strongbc_rhs") - -# Some other useful kernels - -assign_dat_code = """ -void assign_dat(double *dest, double *src) -{ - *dest = *src; -}""" - -assign_dat = op2.Kernel(assign_dat_code, "assign_dat") - -l2norm_diff_sq_code = """ -void l2norm_diff_sq(double *f, double *g, double *norm) -{ - double diff = abs(*f - *g); - *norm += diff*diff; -} -""" - -l2norm_diff_sq = op2.Kernel(l2norm_diff_sq_code, "l2norm_diff_sq") - -# Nonlinear iteration - -# Tol = 1.e-8 -tolsq = 1.e-16 -normsq = op2.Global(1, data=10000.0, name="norm") -solver = op2.Solver() - -while normsq.data[0] > tolsq: - - # Assign result from previous timestep - - op2.par_loop(assign_dat, nodes, - tracer_old(op2.WRITE), - tracer(op2.READ)) - - # Matrix assembly - - mat.zero() - - op2.par_loop(burgers, elements, - mat(op2.INC, (elem_node[op2.i[0]], elem_node[op2.i[1]])), - coords(op2.READ, elem_node), - tracer(op2.READ, elem_node)) - - mat.zero_rows([0, n - 1], 1.0) - - # RHS Assembly - - rhs.zero() - - op2.par_loop(rhs, elements, - b(op2.INC, elem_node[op2.i[0]]), - coords(op2.READ, elem_node), - tracer(op2.READ, elem_node)) - - op2.par_loop(strongbc_rhs, b_nodes, - bdry(op2.READ), - b(op2.WRITE, b_node_node[0])) - - # Solve - - solver.solve(mat, tracer, b) - - # Calculate L2-norm^2 - - normsq = op2.Global(1, data=0.0, name="norm") - op2.par_loop(l2norm_diff_sq, nodes, - tracer(op2.READ), - tracer_old(op2.READ), - normsq(op2.INC)) - - print "L2 Norm squared: %s" % normsq.data[0] - -if opt['plot']: - pylab.plot(coords.data, tracer.data) - pylab.show() From 9a49983f0b187e848309714f6a90a682e5990a71 Mon Sep 17 00:00:00 2001 From: Florian Rathgeber Date: Fri, 4 Apr 2014 15:30:11 +0100 Subject: [PATCH 02/18] Remove stupid MPI demo --- demo/stupid_mpi.py | 172 --------------------------------------------- 1 file changed, 172 deletions(-) delete mode 100644 demo/stupid_mpi.py diff --git a/demo/stupid_mpi.py b/demo/stupid_mpi.py deleted file mode 100644 index 7ac5c9be7..000000000 --- a/demo/stupid_mpi.py +++ /dev/null @@ -1,172 +0,0 @@ -# This file is part of PyOP2 -# -# PyOP2 is Copyright (c) 2012, Imperial College London and -# others. Please see the AUTHORS file in the main source directory for -# a full list of copyright holders. All rights reserved. -# -# Redistribution and use in source and binary forms, with or without -# modification, are permitted provided that the following conditions -# are met: -# -# * Redistributions of source code must retain the above copyright -# notice, this list of conditions and the following disclaimer. -# * Redistributions in binary form must reproduce the above copyright -# notice, this list of conditions and the following disclaimer in the -# documentation and/or other materials provided with the distribution. -# * The name of Imperial College London or that of other -# contributors may not be used to endorse or promote products -# derived from this software without specific prior written -# permission. -# -# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTERS -# ''AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT -# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS -# FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE -# COPYRIGHT HOLDERS OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, -# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES -# (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR -# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) -# HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, -# STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) -# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED -# OF THE POSSIBILITY OF SUCH DAMAGE. - -"""PyOP2 Stupid MPI demo - -This demo repeatedly computes the input mesh geometric center by two means -and scaling the mesh around its center. - -The domain read in from a pickle dump. -""" - -import numpy as np -from numpy.testing import assert_almost_equal, assert_allclose -from cPickle import load -import gzip - -from pyop2 import op2, utils - - -def main(opt): - valuetype = np.float64 - - f = gzip.open(opt['mesh'] + '.' + str(op2.MPI.comm.rank) + '.pickle.gz') - - elements, nodes, elem_node, coords = load(f) - f.close() - coords = op2.Dat(nodes ** 2, coords.data, np.float64, "coords") - varea = op2.Dat(nodes, np.zeros((nodes.total_size, 1), valuetype), valuetype, "varea") - - mesh_center = op2.Kernel("""\ -void -mesh_center(double* coords, double* center, int* count) -{ - center[0] += coords[0]; - center[1] += coords[1]; - *count += 1; -}""", "mesh_center") - - mesh_scale = op2.Kernel("""\ -void -mesh_scale(double* coords, double* center, double* scale) -{ - coords[0] = (coords[0] - center[0]) * scale[0] + center[0]; - coords[1] = (coords[1] - center[1]) * scale[1] + center[1]; -}""", "mesh_scale") - - elem_center = op2.Kernel("""\ -void -elem_center(double* center, double* vcoords[3], int* count) -{ - center[0] += (vcoords[0][0] + vcoords[1][0] + vcoords[2][0]) / 3.0f; - center[1] += (vcoords[0][1] + vcoords[1][1] + vcoords[2][1]) / 3.0f; - *count += 1; -}""", "elem_center") - - dispatch_area = op2.Kernel("""\ -void -dispatch_area(double* vcoords[3], double* area[3]) -{ - double a = 0; - a += vcoords[0][0] * ( vcoords[1][1] - vcoords[2][1] ); - a += vcoords[1][0] * ( vcoords[2][1] - vcoords[0][1] ); - a += vcoords[2][0] * ( vcoords[0][1] - vcoords[1][1] ); - a = fabs(a) / 6.0; - - *area[0] += a; - *area[1] += a; - *area[2] += a; -}""", "dispatch_area") - - collect_area = op2.Kernel("""\ -void -collect_area(double* varea, double* area) -{ - *area += *varea; -}""", "collect_area") - - expected_area = 1.0 - for s in [[1, 2], [2, 1], [3, 3], [2, 5], [5, 2]]: - center1 = op2.Global(2, [0.0, 0.0], valuetype, name='center1') - center2 = op2.Global(2, [0.0, 0.0], valuetype, name='center2') - node_count = op2.Global(1, [0], np.int32, name='node_count') - elem_count = op2.Global(1, [0], np.int32, name='elem_count') - scale = op2.Global(2, s, valuetype, name='scale') - area = op2.Global(1, [0.0], valuetype, name='area') - - op2.par_loop(mesh_center, nodes, - coords(op2.READ), - center1(op2.INC), - node_count(op2.INC)) - center1.data[:] = center1.data[:] / node_count.data[:] - - op2.par_loop(elem_center, elements, - center2(op2.INC), - coords(op2.READ, elem_node), - elem_count(op2.INC)) - center2.data[:] = center2.data[:] / elem_count.data[:] - - op2.par_loop(mesh_scale, nodes, - coords(op2.RW), - center1(op2.READ), - scale(op2.READ)) - - varea.zero() - op2.par_loop(dispatch_area, elements, - coords(op2.READ, elem_node), - varea(op2.INC, elem_node)) - - op2.par_loop(collect_area, nodes, - varea(op2.READ), - area(op2.INC)) - - expected_area *= s[0] * s[1] - - if opt['print_output']: - print "Rank: %d: [%f, %f] [%f, %f] |%f (%f)|" % \ - (op2.MPI.comm.rank, - center1.data[0], center1.data[1], - center2.data[0], center2.data[1], - area.data[0], expected_area) - - if opt['test_output']: - assert_allclose(center1.data, [0.5, 0.5]) - assert_allclose(center2.data, center1.data) - assert_almost_equal(area.data[0], expected_area) - -if __name__ == '__main__': - parser = utils.parser(group=True, description=__doc__) - parser.add_argument('-m', '--mesh', required=True, - help='Base name of mesh pickle \ - (excluding the process number and .pickle extension)') - parser.add_argument('--print-output', action='store_true', help='Print output') - parser.add_argument('--test-output', action='store_true', help='Test output') - - opt = vars(parser.parse_args()) - op2.init(**opt) - - if op2.MPI.comm.size != 3: - print "Stupid demo only works on 3 processes" - op2.MPI.comm.Abort(1) - - main(opt) From d995525c13641c6c9b398ddc1e4200ca72ed2c49 Mon Sep 17 00:00:00 2001 From: Florian Rathgeber Date: Fri, 31 Jan 2014 17:44:47 +0000 Subject: [PATCH 03/18] Remove adv_diff_mpi demo This demo is hard to maintain since it relies on pickled objects that need updating whenever the relevant class changes. Prevent further bitrot by removing it. --- demo/adv_diff_mpi.py | 240 ------------------------------------------- demo/meshes/Makefile | 6 +- 2 files changed, 2 insertions(+), 244 deletions(-) delete mode 100644 demo/adv_diff_mpi.py diff --git a/demo/adv_diff_mpi.py b/demo/adv_diff_mpi.py deleted file mode 100644 index ff9fde529..000000000 --- a/demo/adv_diff_mpi.py +++ /dev/null @@ -1,240 +0,0 @@ -# This file is part of PyOP2 -# -# PyOP2 is Copyright (c) 2012, Imperial College London and -# others. Please see the AUTHORS file in the main source directory for -# a full list of copyright holders. All rights reserved. -# -# Redistribution and use in source and binary forms, with or without -# modification, are permitted provided that the following conditions -# are met: -# -# * Redistributions of source code must retain the above copyright -# notice, this list of conditions and the following disclaimer. -# * Redistributions in binary form must reproduce the above copyright -# notice, this list of conditions and the following disclaimer in the -# documentation and/or other materials provided with the distribution. -# * The name of Imperial College London or that of other -# contributors may not be used to endorse or promote products -# derived from this software without specific prior written -# permission. -# -# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTERS -# ''AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT -# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS -# FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE -# COPYRIGHT HOLDERS OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, -# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES -# (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR -# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) -# HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, -# STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) -# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED -# OF THE POSSIBILITY OF SUCH DAMAGE. - -"""PyOP2 P1 MPI advection-diffusion demo - -This demo solves the advection-diffusion equation by splitting the advection -and diffusion terms. The advection term is advanced in time using an Euler -method and the diffusion term is advanced in time using a theta scheme with -theta = 0.5. - -The domain read in from a pickle dump. - -This demo requires the MAPDES forks of FFC, FIAT and UFL which are found at: - - https://bitbucket.org/mapdes/ffc - https://bitbucket.org/mapdes/fiat - https://bitbucket.org/mapdes/ufl -""" - -import os -import numpy as np -from cPickle import load -import gzip - -from pyop2 import op2, utils -from pyop2.ffc_interface import compile_form -from ufl import * - - -def main(opt): - # Set up finite element problem - - dt = 0.0001 - - T = FiniteElement("Lagrange", "triangle", 1) - V = VectorElement("Lagrange", "triangle", 1) - - p = TrialFunction(T) - q = TestFunction(T) - t = Coefficient(T) - u = Coefficient(V) - a = Coefficient(T) - - diffusivity = 0.1 - - M = p * q * dx - - adv_rhs = (q * t + dt * dot(grad(q), u) * t) * dx - - d = -dt * diffusivity * dot(grad(q), grad(p)) * dx - - diff = M - 0.5 * d - diff_rhs = action(M + 0.5 * d, t) - - # Generate code for mass and rhs assembly. - - adv, = compile_form(M, "adv") - adv_rhs, = compile_form(adv_rhs, "adv_rhs") - diff, = compile_form(diff, "diff") - diff_rhs, = compile_form(diff_rhs, "diff_rhs") - - # Set up simulation data structures - - valuetype = np.float64 - - f = gzip.open(opt['mesh'] + '.' + str(op2.MPI.comm.rank) + '.pickle.gz') - - elements, nodes, elem_node, coords = load(f) - f.close() - coords = op2.Dat(nodes ** 2, coords.data, np.float64, "dcoords") - - num_nodes = nodes.total_size - - sparsity = op2.Sparsity((nodes, nodes), (elem_node, elem_node), "sparsity") - if opt['advection']: - adv_mat = op2.Mat(sparsity, valuetype, "adv_mat") - op2.par_loop(adv, elements, - adv_mat(op2.INC, (elem_node[op2.i[0]], elem_node[op2.i[1]])), - coords(op2.READ, elem_node, flatten=True)) - if opt['diffusion']: - diff_mat = op2.Mat(sparsity, valuetype, "diff_mat") - op2.par_loop(diff, elements, - diff_mat(op2.INC, (elem_node[op2.i[0]], elem_node[op2.i[1]])), - coords(op2.READ, elem_node, flatten=True)) - - tracer_vals = np.zeros(num_nodes, dtype=valuetype) - tracer = op2.Dat(nodes, tracer_vals, valuetype, "tracer") - - b_vals = np.zeros(num_nodes, dtype=valuetype) - b = op2.Dat(nodes, b_vals, valuetype, "b") - - velocity_vals = np.asarray([1.0, 0.0] * num_nodes, dtype=valuetype) - velocity = op2.Dat(nodes ** 2, velocity_vals, valuetype, "velocity") - - # Set initial condition - - i_cond_code = """void i_cond(double *c, double *t) -{ - double A = 0.1; // Normalisation - double D = 0.1; // Diffusivity - double pi = 3.14159265358979; - double x = c[0]-(0.45+%(T)f); - double y = c[1]-0.5; - double r2 = x*x+y*y; - - *t = A*(exp(-r2/(4*D*%(T)f))/(4*pi*D*%(T)f)); -} -""" - - T = 0.01 - - i_cond = op2.Kernel(i_cond_code % {'T': T}, "i_cond") - - op2.par_loop(i_cond, nodes, - coords(op2.READ, flatten=True), - tracer(op2.WRITE)) - - # Assemble and solve - - solver = op2.Solver() - - while T < 0.015: - - # Advection - - if opt['advection']: - b.zero() - op2.par_loop(adv_rhs, elements, - b(op2.INC, elem_node[op2.i[0]]), - coords(op2.READ, elem_node, flatten=True), - tracer(op2.READ, elem_node), - velocity(op2.READ, elem_node)) - - solver.solve(adv_mat, tracer, b) - - # Diffusion - - if opt['diffusion']: - b.zero() - op2.par_loop(diff_rhs, elements, - b(op2.INC, elem_node[op2.i[0]]), - coords(op2.READ, elem_node, flatten=True), - tracer(op2.READ, elem_node)) - - solver.solve(diff_mat, tracer, b) - - T = T + dt - - if opt['print_output'] or opt['test_output']: - analytical_vals = np.zeros(num_nodes, dtype=valuetype) - analytical = op2.Dat(nodes, analytical_vals, valuetype, "analytical") - - i_cond = op2.Kernel(i_cond_code % {'T': T}, "i_cond") - - op2.par_loop(i_cond, nodes, - coords(op2.READ, flatten=True), - analytical(op2.WRITE)) - - # Print error w.r.t. analytical solution - if opt['print_output']: - print "Rank: %d Expected - computed solution: %s" % \ - (op2.MPI.comm.rank, tracer.data - analytical.data) - - if opt['test_output']: - l2norm = dot(t - a, t - a) * dx - l2_kernel, = compile_form(l2norm, "error_norm") - result = op2.Global(1, [0.0]) - op2.par_loop(l2_kernel, elements, - result(op2.INC), - coords(op2.READ, elem_node, flatten=True), - tracer(op2.READ, elem_node), - analytical(op2.READ, elem_node) - ) - if op2.MPI.comm.rank == 0: - with open("adv_diff_mpi.%s.out" % os.path.split(opt['mesh'])[-1], - "w") as out: - out.write(str(result.data[0])) - else: - # hack to prevent mpi communication dangling - result.data - -if __name__ == '__main__': - parser = utils.parser(group=True, description=__doc__) - parser.add_argument('-m', '--mesh', required=True, - help='Base name of mesh pickle \ - (excluding the process number and .pickle extension)') - parser.add_argument('--no-advection', action='store_false', - dest='advection', help='Disable advection') - parser.add_argument('--no-diffusion', action='store_false', - dest='diffusion', help='Disable diffusion') - parser.add_argument('--print-output', action='store_true', help='Print output') - parser.add_argument('-t', '--test-output', action='store_true', - help='Save output for testing') - parser.add_argument('-p', '--profile', action='store_true', - help='Create a cProfile for the run') - - opt = vars(parser.parse_args()) - op2.init(**opt) - - if op2.MPI.comm.size != 3: - print "MPI advection-diffusion demo only works on 3 processes" - op2.MPI.comm.Abort(1) - - if opt['profile']: - import cProfile - filename = 'adv_diff.%s.%d.cprofile' % ( - os.path.split(opt['mesh'])[-1], op2.MPI.comm.rank) - cProfile.run('main(opt)', filename=filename) - else: - main(opt) diff --git a/demo/meshes/Makefile b/demo/meshes/Makefile index bd29e3d91..86ba21b93 100644 --- a/demo/meshes/Makefile +++ b/demo/meshes/Makefile @@ -1,13 +1,11 @@ WGET = wget --no-check-certificate BASEURL = https://spo.doc.ic.ac.uk/meshes/ -PROCS = 0 1 2 -MMS_MESHES = $(foreach mesh, MMS_A MMS_B MMS_C MMS_D, $(foreach proc, $(PROCS), $(mesh).$(proc).pickle.gz)) HDF5_MESHES = new_grid.h5 FE_grid.h5 TRIANGLE_MESHES = $(foreach mesh, small medium large, $(foreach ext, edge ele node, $(mesh).$(ext))) .PHONY : meshes -%.pickle.gz %.h5: +%.h5: $(WGET) $(BASEURL)$@ small.%: @@ -19,4 +17,4 @@ medium.%: large.%: ./generate_mesh large 40 -meshes: $(MMS_MESHES) $(HDF5_MESHES) $(TRIANGLE_MESHES) +meshes: $(HDF5_MESHES) $(TRIANGLE_MESHES) From 3c6aeefa3341354669dd816a18a68549ca9ea6da Mon Sep 17 00:00:00 2001 From: Florian Rathgeber Date: Fri, 4 Apr 2014 15:06:32 +0100 Subject: [PATCH 04/18] Remove regression tests No more need for demo and demo/meshes to be packages. --- .gitignore | 1 - Makefile | 18 ++------ demo/__init__.py | 0 demo/meshes/__init__.py | 0 test/regression/demo | 1 - test/regression/meshes/square.poly | 11 ----- test/regression/test_regression.py | 69 ------------------------------ 7 files changed, 3 insertions(+), 97 deletions(-) delete mode 100644 demo/__init__.py delete mode 100644 demo/meshes/__init__.py delete mode 120000 test/regression/demo delete mode 100644 test/regression/meshes/square.poly delete mode 100644 test/regression/test_regression.py diff --git a/.gitignore b/.gitignore index c6077b025..dced2d737 100644 --- a/.gitignore +++ b/.gitignore @@ -15,4 +15,3 @@ pyop2.log *.node *.geo !cdisk.geo -/test/regression/tests/**/*.out diff --git a/Makefile b/Makefile index 81bf1c431..583dac2de 100644 --- a/Makefile +++ b/Makefile @@ -4,8 +4,6 @@ TEST_BASE_DIR = test UNIT_TEST_DIR = $(TEST_BASE_DIR)/unit -REGRESSION_TEST_DIR = $(TEST_BASE_DIR)/regression - BACKENDS ?= sequential opencl openmp cuda OPENCL_ALL_CTXS := $(shell scripts/detect_opencl_devices) OPENCL_CTXS ?= $(OPENCL_ALL_CTXS) @@ -25,16 +23,14 @@ GIT_REV = $(shell git rev-parse --verify --short HEAD) all: ext -.PHONY : help test lint unit regression doc update_docs ext ext_clean meshes +.PHONY : help test lint unit doc update_docs ext ext_clean meshes help: @echo "make COMMAND with COMMAND one of:" - @echo " test : run lint, unit and regression tests" + @echo " test : run lint and unit tests" @echo " lint : run flake8 code linter" @echo " unit : run unit tests" @echo " unit_BACKEND : run unit tests for BACKEND" - @echo " regression : run regression tests" - @echo " regression_BACKEND : run regression tests for BACKEND" @echo " doc : build sphinx documentation" @echo " serve : launch local web server to serve up documentation" @echo " update_docs : build sphinx documentation and push to GitHub" @@ -44,7 +40,7 @@ help: @echo @echo "Available OpenCL contexts: $(OPENCL_CTXS)" -test: lint unit regression +test: lint unit lint: @flake8 @@ -57,14 +53,6 @@ unit_%: unit_opencl: cd $(UNIT_TEST_DIR); for c in $(OPENCL_CTXS); do PYOPENCL_CTX=$$c $(PYTEST) --backend=opencl; done -regression: $(foreach backend,$(BACKENDS), regression_$(backend)) - -regression_%: - cd $(REGRESSION_TEST_DIR); $(PYTEST) --backend=$* - -regression_opencl: - cd $(REGRESSION_TEST_DIR); for c in $(OPENCL_CTXS); do PYOPENCL_CTX=$$c $(PYTEST) --backend=opencl; done - doc: make -C $(SPHINX_DIR) $(SPHINX_TARGET) SPHINXOPTS=$(SPHINX_OPTS) APIDOCOPTS=$(APIDOC_OPTS) diff --git a/demo/__init__.py b/demo/__init__.py deleted file mode 100644 index e69de29bb..000000000 diff --git a/demo/meshes/__init__.py b/demo/meshes/__init__.py deleted file mode 100644 index e69de29bb..000000000 diff --git a/test/regression/demo b/test/regression/demo deleted file mode 120000 index bf71256cd..000000000 --- a/test/regression/demo +++ /dev/null @@ -1 +0,0 @@ -../../demo \ No newline at end of file diff --git a/test/regression/meshes/square.poly b/test/regression/meshes/square.poly deleted file mode 100644 index b48a8a83c..000000000 --- a/test/regression/meshes/square.poly +++ /dev/null @@ -1,11 +0,0 @@ -4 2 0 0 -1 0 0 -2 1 0 -3 1 1 -4 0 1 -4 1 -1 1 2 3 -2 2 3 2 -3 3 4 3 -4 4 1 1 -0 \ No newline at end of file diff --git a/test/regression/test_regression.py b/test/regression/test_regression.py deleted file mode 100644 index 8fa56487e..000000000 --- a/test/regression/test_regression.py +++ /dev/null @@ -1,69 +0,0 @@ -from os.path import join, dirname, abspath, exists -from subprocess import call - -import numpy as np -import pytest - - -@pytest.fixture(scope='session') -def meshdir(): - return lambda m='': join(join(dirname(abspath(__file__)), 'meshes'), m) - - -@pytest.fixture -def mms_meshes(meshdir): - from demo.meshes.generate_mesh import generate_meshfile - m = [(meshdir('a'), 20), (meshdir('b'), 40), (meshdir('c'), 80), (meshdir('d'), 160)] - for name, layers in m: - if not all(exists(name + ext) for ext in ['.edge', '.ele', '.node']): - generate_meshfile(name, layers) - return m - - -@pytest.fixture -def unstructured_square(meshdir): - m = meshdir('square.1') - if not all(exists(m + ext) for ext in ['.edge', '.ele', '.node']): - call(['triangle', '-e', '-a0.00007717', meshdir('square.poly')]) - return m - - -def test_adv_diff(backend, mms_meshes): - from demo.adv_diff import main, parser - res = np.array([np.sqrt(main(vars(parser.parse_args(['-m', name, '-r'])))) - for name, _ in mms_meshes]) - convergence = np.log2(res[:len(mms_meshes) - 1] / res[1:]) - assert all(convergence > [1.5, 1.85, 1.95]) - - -def test_laplace_ffc(backend): - from demo.laplace_ffc import main, parser - f, x = main(vars(parser.parse_args(['-r']))) - assert sum(abs(f - x)) < 1e-12 - - -def test_mass2d_ffc(backend): - from demo.mass2d_ffc import main, parser - f, x = main(vars(parser.parse_args(['-r']))) - assert sum(abs(f - x)) < 1e-12 - - -def test_mass2d_triangle(backend, unstructured_square): - from demo.mass2d_triangle import main, parser - f, x = main(vars(parser.parse_args(['-m', unstructured_square, '-r']))) - assert np.linalg.norm(f - x) / np.linalg.norm(f) < 1e-6 - - -def test_mass_vector_ffc(backend): - from demo.mass_vector_ffc import main, parser - f, x = main(vars(parser.parse_args(['-r']))) - assert abs(f - x).sum() < 1e-12 - - -@pytest.mark.xfail('config.getvalue("backend")[0] in ("cuda", "opencl")', - reason='Need to expose loops inside conditionals, \ - or to re-design to avoid conditionals') -def test_weak_bcs_ffc(backend): - from demo.weak_bcs_ffc import main, parser - f, x = main(vars(parser.parse_args(['-r']))) - assert abs(f - x).sum() < 1e-12 From 3c096293d2afcb650d757c155566715e8c459eb5 Mon Sep 17 00:00:00 2001 From: Florian Rathgeber Date: Fri, 4 Apr 2014 15:09:55 +0100 Subject: [PATCH 05/18] README: remove mentions of regression tests --- README.rst | 15 ++------------- 1 file changed, 2 insertions(+), 13 deletions(-) diff --git a/README.rst b/README.rst index 1c7cbb6c9..f2e2b33bb 100644 --- a/README.rst +++ b/README.rst @@ -383,19 +383,12 @@ following to ``~/.bashrc`` or ``.env``:: # Add pytest binaries to the path export PATH=${PATH}:${HOME}/.local/bin -The regression tests further require *gmsh* and *triangle*: :: - - sudo apt-get install gmsh triangle-bin unzip - If all tests in our test suite pass, you should be good to go:: make test -This will run both unit and regression tests, the latter require UFL_ and FFC_. - -This will attempt to run tests for all backends and skip those for not -available backends. If the FFC_ fork is not found, tests for the FFC_ interface -are xfailed. +This will run code linting and unit tests, attempting to run for all backends +and skipping those for not available backends. Troubleshooting --------------- @@ -418,10 +411,6 @@ Start with the unit tests with the sequential backend :: py.test test/unit -vsx --tb=short --backend=sequential -Then move on to the regression tests with the sequential backend :: - - py.test test/regression -vsx --tb=short --backend=sequential - With all the sequential tests passing, move on to the next backend in the same manner as required. From b248aafd51602fa9a7322b427e7349dc86f64159 Mon Sep 17 00:00:00 2001 From: Florian Rathgeber Date: Fri, 4 Apr 2014 14:30:49 +0100 Subject: [PATCH 06/18] Remove advection-diffusion demo --- demo/adv_diff.py | 246 ----------------------------------------------- 1 file changed, 246 deletions(-) delete mode 100644 demo/adv_diff.py diff --git a/demo/adv_diff.py b/demo/adv_diff.py deleted file mode 100644 index 511839e8f..000000000 --- a/demo/adv_diff.py +++ /dev/null @@ -1,246 +0,0 @@ -# This file is part of PyOP2 -# -# PyOP2 is Copyright (c) 2012, Imperial College London and -# others. Please see the AUTHORS file in the main source directory for -# a full list of copyright holders. All rights reserved. -# -# Redistribution and use in source and binary forms, with or without -# modification, are permitted provided that the following conditions -# are met: -# -# * Redistributions of source code must retain the above copyright -# notice, this list of conditions and the following disclaimer. -# * Redistributions in binary form must reproduce the above copyright -# notice, this list of conditions and the following disclaimer in the -# documentation and/or other materials provided with the distribution. -# * The name of Imperial College London or that of other -# contributors may not be used to endorse or promote products -# derived from this software without specific prior written -# permission. -# -# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTERS -# ''AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT -# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS -# FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE -# COPYRIGHT HOLDERS OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, -# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES -# (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR -# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) -# HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, -# STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) -# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED -# OF THE POSSIBILITY OF SUCH DAMAGE. - -"""PyOP2 P1 advection-diffusion with operator splitting demo - -This demo solves the advection-diffusion equation by splitting the advection -and diffusion terms. The advection term is advanced in time using an Euler -method and the diffusion term is advanced in time using a theta scheme with -theta = 0.5. - -The domain read in from a triangle file. - -This demo requires the MAPDES forks of FFC, FIAT and UFL which are found at: - - https://bitbucket.org/mapdes/ffc - https://bitbucket.org/mapdes/fiat - https://bitbucket.org/mapdes/ufl - -FEniCS Viper is optionally used to visualise the solution. -""" -import os -import numpy as np - -from pyop2 import op2, utils -from pyop2.ffc_interface import compile_form -from triangle_reader import read_triangle -from ufl import * - - -def viper_shape(array): - """Flatten a numpy array into one dimension to make it suitable for - passing to Viper.""" - return array.reshape((array.shape[0])) - - -def main(opt): - # Set up finite element problem - - dt = 0.0001 - - T = FiniteElement("Lagrange", "triangle", 1) - V = VectorElement("Lagrange", "triangle", 1) - - p = TrialFunction(T) - q = TestFunction(T) - t = Coefficient(T) - u = Coefficient(V) - a = Coefficient(T) - - diffusivity = 0.1 - - M = p * q * dx - - adv_rhs = (q * t + dt * dot(grad(q), u) * t) * dx - - d = -dt * diffusivity * dot(grad(q), grad(p)) * dx - - diff = M - 0.5 * d - diff_rhs = action(M + 0.5 * d, t) - - # Generate code for mass and rhs assembly. - - adv, = compile_form(M, "adv") - adv_rhs, = compile_form(adv_rhs, "adv_rhs") - diff, = compile_form(diff, "diff") - diff_rhs, = compile_form(diff_rhs, "diff_rhs") - - # Set up simulation data structures - - valuetype = np.float64 - - nodes, coords, elements, elem_node = read_triangle(opt['mesh']) - - num_nodes = nodes.size - - sparsity = op2.Sparsity((nodes, nodes), (elem_node, elem_node), "sparsity") - if opt['advection']: - adv_mat = op2.Mat(sparsity, valuetype, "adv_mat") - op2.par_loop(adv, elements, - adv_mat(op2.INC, (elem_node[op2.i[0]], elem_node[op2.i[1]])), - coords(op2.READ, elem_node, flatten=True)) - if opt['diffusion']: - diff_mat = op2.Mat(sparsity, valuetype, "diff_mat") - op2.par_loop(diff, elements, - diff_mat(op2.INC, (elem_node[op2.i[0]], elem_node[op2.i[1]])), - coords(op2.READ, elem_node, flatten=True)) - - tracer_vals = np.zeros(num_nodes, dtype=valuetype) - tracer = op2.Dat(nodes, tracer_vals, valuetype, "tracer") - - b_vals = np.zeros(num_nodes, dtype=valuetype) - b = op2.Dat(nodes, b_vals, valuetype, "b") - - velocity_vals = np.asarray([1.0, 0.0] * num_nodes, dtype=valuetype) - velocity = op2.Dat(nodes ** 2, velocity_vals, valuetype, "velocity") - - # Set initial condition - - i_cond_code = """void i_cond(double *c, double *t) -{ - double A = 0.1; // Normalisation - double D = 0.1; // Diffusivity - double pi = 3.14159265358979; - double x = c[0]-(0.45+%(T)f); - double y = c[1]-0.5; - double r2 = x*x+y*y; - - *t = A*(exp(-r2/(4*D*%(T)f))/(4*pi*D*%(T)f)); -} -""" - - T = 0.01 - - i_cond = op2.Kernel(i_cond_code % {'T': T}, "i_cond") - - op2.par_loop(i_cond, nodes, - coords(op2.READ, flatten=True), - tracer(op2.WRITE)) - - # Assemble and solve - if opt['visualize']: - vis_coords = np.asarray([[x, y, 0.0] for x, y in coords.data_ro], - dtype=np.float64) - import viper - v = viper.Viper(x=viper_shape(tracer.data_ro), - coordinates=vis_coords, cells=elem_node.values) - - solver = op2.Solver() - - while T < 0.015: - - # Advection - - if opt['advection']: - b.zero() - op2.par_loop(adv_rhs, elements, - b(op2.INC, elem_node[op2.i[0]]), - coords(op2.READ, elem_node, flatten=True), - tracer(op2.READ, elem_node), - velocity(op2.READ, elem_node, flatten=True)) - - solver.solve(adv_mat, tracer, b) - - # Diffusion - - if opt['diffusion']: - b.zero() - op2.par_loop(diff_rhs, elements, - b(op2.INC, elem_node[op2.i[0]]), - coords(op2.READ, elem_node, flatten=True), - tracer(op2.READ, elem_node)) - - solver.solve(diff_mat, tracer, b) - - if opt['visualize']: - v.update(viper_shape(tracer.data_ro)) - - T = T + dt - - if opt['print_output'] or opt['test_output'] or opt['return_output']: - analytical_vals = np.zeros(num_nodes, dtype=valuetype) - analytical = op2.Dat(nodes, analytical_vals, valuetype, "analytical") - - i_cond = op2.Kernel(i_cond_code % {'T': T}, "i_cond") - - op2.par_loop(i_cond, nodes, - coords(op2.READ, flatten=True), - analytical(op2.WRITE)) - - # Print error w.r.t. analytical solution - if opt['print_output']: - print "Expected - computed solution: %s" % (tracer.data - analytical.data) - - if opt['test_output'] or opt['return_output']: - l2norm = dot(t - a, t - a) * dx - l2_kernel, = compile_form(l2norm, "error_norm") - result = op2.Global(1, [0.0]) - op2.par_loop(l2_kernel, elements, - result(op2.INC), - coords(op2.READ, elem_node, flatten=True), - tracer(op2.READ, elem_node), - analytical(op2.READ, elem_node)) - if opt['test_output']: - with open("adv_diff.%s.out" % os.path.split(opt['mesh'])[-1], "w") as out: - out.write(str(result.data[0]) + "\n") - if opt['return_output']: - return result.data[0] - -parser = utils.parser(group=True, description=__doc__) -parser.add_argument('-m', '--mesh', required=True, - help='Base name of triangle mesh \ - (excluding the .ele or .node extension)') -parser.add_argument('-v', '--visualize', action='store_true', - help='Visualize the result using viper') -parser.add_argument('--no-advection', action='store_false', - dest='advection', help='Disable advection') -parser.add_argument('--no-diffusion', action='store_false', - dest='diffusion', help='Disable diffusion') -parser.add_argument('--print-output', action='store_true', help='Print output') -parser.add_argument('-r', '--return-output', action='store_true', - help='Return output for testing') -parser.add_argument('-t', '--test-output', action='store_true', - help='Save output for testing') -parser.add_argument('-p', '--profile', action='store_true', - help='Create a cProfile for the run') - -if __name__ == '__main__': - opt = vars(parser.parse_args()) - op2.init(**opt) - - if opt['profile']: - import cProfile - filename = 'adv_diff.%s.cprofile' % os.path.split(opt['mesh'])[-1] - cProfile.run('main(opt)', filename=filename) - else: - main(opt) From a23a4e27bede36e05a47279d65925ac6f5b69af4 Mon Sep 17 00:00:00 2001 From: Florian Rathgeber Date: Fri, 4 Apr 2014 14:31:14 +0100 Subject: [PATCH 07/18] Remove non-split advection-diffusion demo --- demo/adv_diff_nonsplit.py | 177 -------------------------------------- 1 file changed, 177 deletions(-) delete mode 100644 demo/adv_diff_nonsplit.py diff --git a/demo/adv_diff_nonsplit.py b/demo/adv_diff_nonsplit.py deleted file mode 100644 index 1900d3a6e..000000000 --- a/demo/adv_diff_nonsplit.py +++ /dev/null @@ -1,177 +0,0 @@ -# This file is part of PyOP2 -# -# PyOP2 is Copyright (c) 2012, Imperial College London and -# others. Please see the AUTHORS file in the main source directory for -# a full list of copyright holders. All rights reserved. -# -# Redistribution and use in source and binary forms, with or without -# modification, are permitted provided that the following conditions -# are met: -# -# * Redistributions of source code must retain the above copyright -# notice, this list of conditions and the following disclaimer. -# * Redistributions in binary form must reproduce the above copyright -# notice, this list of conditions and the following disclaimer in the -# documentation and/or other materials provided with the distribution. -# * The name of Imperial College London or that of other -# contributors may not be used to endorse or promote products -# derived from this software without specific prior written -# permission. -# -# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTERS -# ''AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT -# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS -# FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE -# COPYRIGHT HOLDERS OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, -# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES -# (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR -# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) -# HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, -# STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) -# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED -# OF THE POSSIBILITY OF SUCH DAMAGE. - -"""PyOP2 P1 advection-diffusion demo - -This demo solves the advection-diffusion equation and is advanced in time using -a theta scheme with theta = 0.5. - -The domain read in from a triangle file. - -This demo requires the MAPDES forks of FFC, FIAT and UFL which are found at: - - https://bitbucket.org/mapdes/ffc - https://bitbucket.org/mapdes/fiat - https://bitbucket.org/mapdes/ufl - -FEniCS Viper is optionally used to visualise the solution. -""" - -from pyop2 import op2, utils -from pyop2.ffc_interface import compile_form -from triangle_reader import read_triangle -from ufl import * - -import numpy as np - - -def viper_shape(array): - """Flatten a numpy array into one dimension to make it suitable for - passing to Viper.""" - return array.reshape((array.shape[0])) - -parser = utils.parser(group=True, description=__doc__) -parser.add_argument('-m', '--mesh', required=True, - help='Base name of triangle mesh \ - (excluding the .ele or .node extension)') -parser.add_argument('-v', '--visualize', action='store_true', - help='Visualize the result using viper') -opt = vars(parser.parse_args()) -op2.init(**opt) - -# Set up finite element problem - -dt = 0.0001 - -T = FiniteElement("Lagrange", "triangle", 1) -V = VectorElement("Lagrange", "triangle", 1) - -p = TrialFunction(T) -q = TestFunction(T) -t = Coefficient(T) -u = Coefficient(V) - -diffusivity = 0.1 - -M = p * q * dx - -d = dt * (diffusivity * dot(grad(q), grad(p)) - dot(grad(q), u) * p) * dx - -a = M + 0.5 * d -L = action(M - 0.5 * d, t) - -# Generate code for mass and rhs assembly. - -lhs, = compile_form(a, "lhs") -rhs, = compile_form(L, "rhs") - -# Set up simulation data structures - -valuetype = np.float64 - -nodes, coords, elements, elem_node = read_triangle(opt['mesh']) - -num_nodes = nodes.size - -sparsity = op2.Sparsity((nodes, nodes), (elem_node, elem_node), "sparsity") -mat = op2.Mat(sparsity, valuetype, "mat") - -tracer_vals = np.zeros(num_nodes, dtype=valuetype) -tracer = op2.Dat(nodes, tracer_vals, valuetype, "tracer") - -b_vals = np.zeros(num_nodes, dtype=valuetype) -b = op2.Dat(nodes, b_vals, valuetype, "b") - -velocity_vals = np.asarray([1.0, 0.0] * num_nodes, dtype=valuetype) -velocity = op2.Dat(nodes ** 2, velocity_vals, valuetype, "velocity") - -# Set initial condition - -i_cond_code = """ -void i_cond(double *c, double *t) -{ - double i_t = 0.1; // Initial time - double A = 0.1; // Normalisation - double D = 0.1; // Diffusivity - double pi = 3.141459265358979; - double x = c[0]-0.5; - double y = c[1]-0.5; - double r = sqrt(x*x+y*y); - - if (r<0.25) - *t = A*(exp((-(r*r))/(4*D*i_t))/(4*pi*D*i_t)); - else - *t = 0.0; -} -""" - -i_cond = op2.Kernel(i_cond_code, "i_cond") - -op2.par_loop(i_cond, nodes, - coords(op2.READ, flatten=True), - tracer(op2.WRITE)) - -# Assemble and solve - -T = 0.1 - -if opt['visualize']: - vis_coords = np.asarray([[x, y, 0.0] for x, y in coords.data_ro], - dtype=np.float64) - import viper - v = viper.Viper(x=viper_shape(tracer.data_ro), - coordinates=vis_coords, cells=elem_node.values) - -solver = op2.Solver() - -while T < 0.2: - - mat.zero() - op2.par_loop(lhs, elements, - mat(op2.INC, (elem_node[op2.i[0]], elem_node[op2.i[1]])), - coords(op2.READ, elem_node, flatten=True), - velocity(op2.READ, elem_node)) - - b.zero() - op2.par_loop(rhs, elements, - b(op2.INC, elem_node[op2.i[0]]), - coords(op2.READ, elem_node, flatten=True), - tracer(op2.READ, elem_node), - velocity(op2.READ, elem_node)) - - solver.solve(mat, tracer, b) - - if opt['visualize']: - v.update(viper_shape(tracer.data_ro)) - - T = T + dt From ad73e0ad48f307687c4ef2b593cee6810c5e6f61 Mon Sep 17 00:00:00 2001 From: Florian Rathgeber Date: Fri, 4 Apr 2014 14:31:42 +0100 Subject: [PATCH 08/18] Remove laplace demo --- demo/laplace_ffc.py | 172 -------------------------------------------- 1 file changed, 172 deletions(-) delete mode 100644 demo/laplace_ffc.py diff --git a/demo/laplace_ffc.py b/demo/laplace_ffc.py deleted file mode 100644 index 88c1c73e2..000000000 --- a/demo/laplace_ffc.py +++ /dev/null @@ -1,172 +0,0 @@ -# This file is part of PyOP2 -# -# PyOP2 is Copyright (c) 2012, Imperial College London and -# others. Please see the AUTHORS file in the main source directory for -# a full list of copyright holders. All rights reserved. -# -# Redistribution and use in source and binary forms, with or without -# modification, are permitted provided that the following conditions -# are met: -# -# * Redistributions of source code must retain the above copyright -# notice, this list of conditions and the following disclaimer. -# * Redistributions in binary form must reproduce the above copyright -# notice, this list of conditions and the following disclaimer in the -# documentation and/or other materials provided with the distribution. -# * The name of Imperial College London or that of other -# contributors may not be used to endorse or promote products -# derived from this software without specific prior written -# permission. -# -# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTERS -# ''AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT -# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS -# FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE -# COPYRIGHT HOLDERS OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, -# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES -# (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR -# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) -# HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, -# STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) -# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED -# OF THE POSSIBILITY OF SUCH DAMAGE. - -"""PyOP2 laplace equation demo - -This demo uses ffc-generated kernels to solve the Laplace equation on a unit -square with boundary conditions: - - u = 1 on y = 0 - u = 2 on y = 1 - -The domain is meshed as follows: - - *-*-* - |/|/| - *-*-* - |/|/| - *-*-* - -This demo requires the MAPDES forks of FFC, FIAT and UFL which are found at: - - https://bitbucket.org/mapdes/ffc - https://bitbucket.org/mapdes/fiat - https://bitbucket.org/mapdes/ufl -""" - -from pyop2 import op2, utils -from pyop2.ffc_interface import compile_form -from ufl import * - -import numpy as np - - -def main(opt): - # Set up finite element problem - - E = FiniteElement("Lagrange", "triangle", 1) - - v = TestFunction(E) - u = TrialFunction(E) - f = Coefficient(E) - - a = dot(grad(v,), grad(u)) * dx - L = v * f * dx - - # Generate code for Laplacian and rhs assembly. - - laplacian, = compile_form(a, "laplacian") - rhs, = compile_form(L, "rhs") - - # Set up simulation data structures - - NUM_ELE = 8 - NUM_NODES = 9 - NUM_BDRY_NODE = 6 - valuetype = np.float64 - - nodes = op2.Set(NUM_NODES, "nodes") - elements = op2.Set(NUM_ELE, "elements") - bdry_nodes = op2.Set(NUM_BDRY_NODE, "boundary_nodes") - - elem_node_map = np.array([0, 1, 4, 4, 3, 0, 1, 2, 5, 5, 4, 1, 3, 4, 7, 7, - 6, 3, 4, 5, 8, 8, 7, 4], dtype=np.uint32) - elem_node = op2.Map(elements, nodes, 3, elem_node_map, "elem_node") - - bdry_node_node_map = np.array([0, 1, 2, 6, 7, 8], dtype=valuetype) - bdry_node_node = op2.Map(bdry_nodes, nodes, 1, bdry_node_node_map, - "bdry_node_node") - - sparsity = op2.Sparsity((nodes, nodes), (elem_node, elem_node), "sparsity") - mat = op2.Mat(sparsity, valuetype, "mat") - - coord_vals = np.array([(0.0, 0.0), (0.5, 0.0), (1.0, 0.0), - (0.0, 0.5), (0.5, 0.5), (1.0, 0.5), - (0.0, 1.0), (0.5, 1.0), (1.0, 1.0)], - dtype=valuetype) - coords = op2.Dat(nodes ** 2, coord_vals, valuetype, "coords") - - u_vals = np.array([1.0, 1.0, 1.0, 1.5, 1.5, 1.5, 2.0, 2.0, 2.0]) - f = op2.Dat(nodes, np.zeros(NUM_NODES, dtype=valuetype), valuetype, "f") - b = op2.Dat(nodes, np.zeros(NUM_NODES, dtype=valuetype), valuetype, "b") - x = op2.Dat(nodes, np.zeros(NUM_NODES, dtype=valuetype), valuetype, "x") - u = op2.Dat(nodes, u_vals, valuetype, "u") - - bdry_vals = np.array([1.0, 1.0, 1.0, 2.0, 2.0, 2.0], dtype=valuetype) - bdry = op2.Dat(bdry_nodes, bdry_vals, valuetype, "bdry") - - # Assemble matrix and rhs - - op2.par_loop(laplacian, elements, - mat(op2.INC, (elem_node[op2.i[0]], elem_node[op2.i[1]])), - coords(op2.READ, elem_node, flatten=True)) - - op2.par_loop(rhs, elements, - b(op2.INC, elem_node[op2.i[0]]), - coords(op2.READ, elem_node, flatten=True), - f(op2.READ, elem_node)) - - # Apply strong BCs - - mat.zero_rows([0, 1, 2, 6, 7, 8], 1.0) - strongbc_rhs = op2.Kernel(""" - void strongbc_rhs(double *val, double *target) { *target = *val; } - """, "strongbc_rhs") - op2.par_loop(strongbc_rhs, bdry_nodes, - bdry(op2.READ), - b(op2.WRITE, bdry_node_node[0])) - - solver = op2.Solver(ksp_type='gmres') - solver.solve(mat, x, b) - - # Print solution - if opt['print_output']: - print "Expected solution: %s" % u.data - print "Computed solution: %s" % x.data - - # Save output (if necessary) - if opt['return_output']: - return u.data, x.data - if opt['save_output']: - import pickle - with open("laplace.out", "w") as out: - pickle.dump((u.data, x.data), out) - -parser = utils.parser(group=True, description=__doc__) -parser.add_argument('--print-output', action='store_true', help='Print output') -parser.add_argument('-r', '--return-output', action='store_true', - help='Return output for testing') -parser.add_argument('-s', '--save-output', action='store_true', - help='Save output for testing') -parser.add_argument('-p', '--profile', action='store_true', - help='Create a cProfile for the run') - -if __name__ == '__main__': - opt = vars(parser.parse_args()) - op2.init(**opt) - - if opt['profile']: - import cProfile - cProfile.run('main(opt)', filename='laplace_ffc.cprofile') - else: - main(opt) From be7c681910c043d2f786c270f4bdcbca975d374f Mon Sep 17 00:00:00 2001 From: Florian Rathgeber Date: Fri, 4 Apr 2014 14:32:04 +0100 Subject: [PATCH 09/18] Remove 2D mass demo --- demo/mass2d_ffc.py | 137 --------------------------------------------- 1 file changed, 137 deletions(-) delete mode 100644 demo/mass2d_ffc.py diff --git a/demo/mass2d_ffc.py b/demo/mass2d_ffc.py deleted file mode 100644 index 18335023f..000000000 --- a/demo/mass2d_ffc.py +++ /dev/null @@ -1,137 +0,0 @@ -# This file is part of PyOP2 -# -# PyOP2 is Copyright (c) 2012, Imperial College London and -# others. Please see the AUTHORS file in the main source directory for -# a full list of copyright holders. All rights reserved. -# -# Redistribution and use in source and binary forms, with or without -# modification, are permitted provided that the following conditions -# are met: -# -# * Redistributions of source code must retain the above copyright -# notice, this list of conditions and the following disclaimer. -# * Redistributions in binary form must reproduce the above copyright -# notice, this list of conditions and the following disclaimer in the -# documentation and/or other materials provided with the distribution. -# * The name of Imperial College London or that of other -# contributors may not be used to endorse or promote products -# derived from this software without specific prior written -# permission. -# -# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTERS -# ''AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT -# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS -# FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE -# COPYRIGHT HOLDERS OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, -# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES -# (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR -# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) -# HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, -# STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) -# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED -# OF THE POSSIBILITY OF SUCH DAMAGE. - -"""PyOP2 2D mass equation demo - -This is a demo of the use of ffc to generate kernels. It solves the identity -equation on a quadrilateral domain. - -This demo requires the MAPDES forks of FFC, FIAT and UFL which are found at: - - https://bitbucket.org/mapdes/ffc - https://bitbucket.org/mapdes/fiat - https://bitbucket.org/mapdes/ufl -""" - -from pyop2 import op2, utils -from pyop2.ffc_interface import compile_form -from ufl import * -import numpy as np - - -def main(opt): - # Set up finite element identity problem - - E = FiniteElement("Lagrange", "triangle", 1) - - v = TestFunction(E) - u = TrialFunction(E) - f = Coefficient(E) - - a = v * u * dx - L = v * f * dx - - # Generate code for mass and rhs assembly. - - mass, = compile_form(a, "mass") - rhs, = compile_form(L, "rhs") - - # Set up simulation data structures - - NUM_ELE = 2 - NUM_NODES = 4 - valuetype = np.float64 - - nodes = op2.Set(NUM_NODES, "nodes") - elements = op2.Set(NUM_ELE, "elements") - - elem_node_map = np.array([0, 1, 3, 2, 3, 1], dtype=np.uint32) - elem_node = op2.Map(elements, nodes, 3, elem_node_map, "elem_node") - - sparsity = op2.Sparsity((nodes, nodes), (elem_node, elem_node), "sparsity") - mat = op2.Mat(sparsity, valuetype, "mat") - - coord_vals = np.array([(0.0, 0.0), (2.0, 0.0), (1.0, 1.0), (0.0, 1.5)], - dtype=valuetype) - coords = op2.Dat(nodes ** 2, coord_vals, valuetype, "coords") - - f = op2.Dat(nodes, np.array([1.0, 2.0, 3.0, 4.0]), valuetype, "f") - b = op2.Dat(nodes, np.zeros(NUM_NODES, dtype=valuetype), valuetype, "b") - x = op2.Dat(nodes, np.zeros(NUM_NODES, dtype=valuetype), valuetype, "x") - - # Assemble and solve - - op2.par_loop(mass, elements, - mat(op2.INC, (elem_node[op2.i[0]], elem_node[op2.i[1]])), - coords(op2.READ, elem_node, flatten=True)) - - op2.par_loop(rhs, elements, - b(op2.INC, elem_node[op2.i[0]]), - coords(op2.READ, elem_node, flatten=True), - f(op2.READ, elem_node)) - - solver = op2.Solver() - solver.solve(mat, x, b) - - # Print solution - if opt['print_output']: - print "Expected solution: %s" % f.data - print "Computed solution: %s" % x.data - - # Save output (if necessary) - if opt['return_output']: - return f.data, x.data - if opt['save_output']: - import pickle - with open("mass2d.out", "w") as out: - pickle.dump((f.data, x.data), out) - -parser = utils.parser(group=True, description=__doc__) -parser.add_argument('--print-output', action='store_true', help='Print output') -parser.add_argument('-r', '--return-output', action='store_true', - help='Return output for testing') -parser.add_argument('-s', '--save-output', - action='store_true', - help='Save the output of the run (used for testing)') -parser.add_argument('-p', '--profile', action='store_true', - help='Create a cProfile for the run') - -if __name__ == '__main__': - opt = vars(parser.parse_args()) - op2.init(**opt) - - if opt['profile']: - import cProfile - cProfile.run('main(opt)', filename='mass2d_ffc.cprofile') - else: - main(opt) From 4fd078a6d717fa366d3955670c49402be2ddbc24 Mon Sep 17 00:00:00 2001 From: Florian Rathgeber Date: Fri, 4 Apr 2014 14:32:32 +0100 Subject: [PATCH 10/18] Remove 2D mass MPI demo --- demo/mass2d_mpi.py | 166 --------------------------------------------- 1 file changed, 166 deletions(-) delete mode 100644 demo/mass2d_mpi.py diff --git a/demo/mass2d_mpi.py b/demo/mass2d_mpi.py deleted file mode 100644 index 90bec6271..000000000 --- a/demo/mass2d_mpi.py +++ /dev/null @@ -1,166 +0,0 @@ -# This file is part of PyOP2 -# -# PyOP2 is Copyright (c) 2012, Imperial College London and -# others. Please see the AUTHORS file in the main source directory for -# a full list of copyright holders. All rights reserved. -# -# Redistribution and use in source and binary forms, with or without -# modification, are permitted provided that the following conditions -# are met: -# -# * Redistributions of source code must retain the above copyright -# notice, this list of conditions and the following disclaimer. -# * Redistributions in binary form must reproduce the above copyright -# notice, this list of conditions and the following disclaimer in the -# documentation and/or other materials provided with the distribution. -# * The name of Imperial College London or that of other -# contributors may not be used to endorse or promote products -# derived from this software without specific prior written -# permission. -# -# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTERS -# ''AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT -# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS -# FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE -# COPYRIGHT HOLDERS OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, -# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES -# (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR -# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) -# HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, -# STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) -# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED -# OF THE POSSIBILITY OF SUCH DAMAGE. - -"""PyOP2 2D mass equation demo (MPI version) - -This is a demo of the use of ffc to generate kernels. It solves the identity -equation on a quadrilateral domain. - -This demo requires the MAPDES forks of FFC, FIAT and UFL which are found at: - - https://bitbucket.org/mapdes/ffc - https://bitbucket.org/mapdes/fiat - https://bitbucket.org/mapdes/ufl -""" - -from pyop2 import op2, utils -from pyop2.ffc_interface import compile_form -from ufl import * -import numpy as np -from petsc4py import PETSc - -parser = utils.parser(group=True, description=__doc__) -parser.add_argument('-s', '--save-output', - action='store_true', - help='Save the output of the run') -parser.add_argument('-t', '--test-output', - action='store_true', - help='Save output for testing') -opt = vars(parser.parse_args()) -op2.init(**opt) - -# Set up finite element identity problem - -E = FiniteElement("Lagrange", "triangle", 1) - -v = TestFunction(E) -u = TrialFunction(E) -f = Coefficient(E) - -a = v * u * dx -L = v * f * dx - -# Generate code for mass and rhs assembly. - -mass, = compile_form(a, "mass") -rhs, = compile_form(L, "rhs") - -# Set up simulation data structures - -NUM_ELE = (0, 1, 2, 2) -NUM_NODES = (0, 2, 4, 4) -valuetype = np.float64 - -if op2.MPI.comm.size != 2: - print "MPI mass2d demo only works on two processes" - op2.MPI.comm.Abort(1) - -if op2.MPI.comm.rank == 0: - node_global_to_universal = np.asarray([0, 1, 2, 3], dtype=PETSc.IntType) - node_halo = op2.Halo(sends={1: [0, 1]}, receives={1: [2, 3]}, - gnn2unn=node_global_to_universal) - element_halo = op2.Halo(sends={1: [0]}, receives={1: [1]}) -elif op2.MPI.comm.rank == 1: - node_global_to_universal = np.asarray([2, 3, 1, 0], dtype=PETSc.IntType) - node_halo = op2.Halo(sends={0: [0, 1]}, receives={0: [3, 2]}, - gnn2unn=node_global_to_universal) - element_halo = op2.Halo(sends={0: [0]}, receives={0: [1]}) -else: - op2.MPI.comm.Abort(1) -nodes = op2.Set(NUM_NODES, "nodes", halo=node_halo) -elements = op2.Set(NUM_ELE, "elements", halo=element_halo) - -if op2.MPI.comm.rank == 0: - elem_node_map = np.asarray([0, 1, 3, 2, 3, 1], dtype=np.uint32) -elif op2.MPI.comm.rank == 1: - elem_node_map = np.asarray([0, 1, 2, 2, 3, 1], dtype=np.uint32) -else: - op2.MPI.comm.Abort(1) - -elem_node = op2.Map(elements, nodes, 3, elem_node_map, "elem_node") - -sparsity = op2.Sparsity((nodes, nodes), (elem_node, elem_node), "sparsity") -mat = op2.Mat(sparsity, valuetype, "mat") - -if op2.MPI.comm.rank == 0: - coord_vals = np.asarray([(0.0, 0.0), (2.0, 0.0), (1.0, 1.0), (0.0, 1.5)], - dtype=valuetype) -elif op2.MPI.comm.rank == 1: - coord_vals = np.asarray([(1, 1), (0, 1.5), (2, 0), (0, 0)], - dtype=valuetype) -else: - op2.MPI.comm.Abort(1) -coords = op2.Dat(nodes ** 2, coord_vals, valuetype, "coords") - -if op2.MPI.comm.rank == 0: - f_vals = np.asarray([1.0, 2.0, 3.0, 4.0], dtype=valuetype) -elif op2.MPI.comm.rank == 1: - f_vals = np.asarray([3.0, 4.0, 2.0, 1.0], dtype=valuetype) -else: - op2.MPI.comm.Abort(1) -b_vals = np.asarray([0.0] * NUM_NODES[3], dtype=valuetype) -x_vals = np.asarray([0.0] * NUM_NODES[3], dtype=valuetype) -f = op2.Dat(nodes, f_vals, valuetype, "f") -b = op2.Dat(nodes, b_vals, valuetype, "b") -x = op2.Dat(nodes, x_vals, valuetype, "x") - -# Assemble and solve - -op2.par_loop(mass, elements, - mat(op2.INC, (elem_node[op2.i[0]], elem_node[op2.i[1]])), - coords(op2.READ, elem_node, flatten=True)) - -op2.par_loop(rhs, elements, - b(op2.INC, elem_node[op2.i[0]]), - coords(op2.READ, elem_node, flatten=True), - f(op2.READ, elem_node)) - -solver = op2.Solver() -solver.solve(mat, x, b) - - -# Compute error in solution -error = (f.data[:f.dataset.size] - x.data[:x.dataset.size]) - -# Print error solution -print "Rank: %d Expected - computed solution: %s" % \ - (op2.MPI.comm.rank, error) - -# Save output (if necessary) -if opt['save_output']: - raise RuntimeException('Writing distributed Dats not yet supported') - -if opt['test_output']: - import pickle - with open("mass2d_mpi_%d.out" % op2.MPI.comm.rank, "w") as out: - pickle.dump(error, out) From adff6a0cd0a2565eab7b413820e1c2868ba2cae6 Mon Sep 17 00:00:00 2001 From: Florian Rathgeber Date: Fri, 4 Apr 2014 14:33:02 +0100 Subject: [PATCH 11/18] Remove 2D mass triangle demo --- demo/mass2d_triangle.py | 139 ---------------------------------------- 1 file changed, 139 deletions(-) delete mode 100644 demo/mass2d_triangle.py diff --git a/demo/mass2d_triangle.py b/demo/mass2d_triangle.py deleted file mode 100644 index 9dcd7b544..000000000 --- a/demo/mass2d_triangle.py +++ /dev/null @@ -1,139 +0,0 @@ -# This file is part of PyOP2 -# -# PyOP2 is Copyright (c) 2012, Imperial College London and -# others. Please see the AUTHORS file in the main source directory for -# a full list of copyright holders. All rights reserved. -# -# Redistribution and use in source and binary forms, with or without -# modification, are permitted provided that the following conditions -# are met: -# -# * Redistributions of source code must retain the above copyright -# notice, this list of conditions and the following disclaimer. -# * Redistributions in binary form must reproduce the above copyright -# notice, this list of conditions and the following disclaimer in the -# documentation and/or other materials provided with the distribution. -# * The name of Imperial College London or that of other -# contributors may not be used to endorse or promote products -# derived from this software without specific prior written -# permission. -# -# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTERS -# ''AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT -# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS -# FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE -# COPYRIGHT HOLDERS OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, -# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES -# (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR -# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) -# HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, -# STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) -# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED -# OF THE POSSIBILITY OF SUCH DAMAGE. - -"""PyOP2 2D mass equation demo - -This demo solves the identity equation on a domain read in from a triangle -file. - -This demo requires the MAPDES forks of FFC, FIAT and UFL which are found at: - - https://bitbucket.org/mapdes/ffc - https://bitbucket.org/mapdes/fiat - https://bitbucket.org/mapdes/ufl -""" - -from pyop2 import op2, utils -from pyop2.ffc_interface import compile_form -from triangle_reader import read_triangle -from ufl import * - -import numpy as np - - -def main(opt): - # Set up finite element identity problem - - E = FiniteElement("Lagrange", "triangle", 1) - - v = TestFunction(E) - u = TrialFunction(E) - f = Coefficient(E) - - a = v * u * dx - L = v * f * dx - - # Generate code for mass and rhs assembly. - - mass, = compile_form(a, "mass") - rhs, = compile_form(L, "rhs") - - # Set up simulation data structures - - valuetype = np.float64 - - nodes, coords, elements, elem_node = read_triangle(opt['mesh']) - - sparsity = op2.Sparsity((nodes, nodes), (elem_node, elem_node), "sparsity") - mat = op2.Mat(sparsity, valuetype, "mat") - - b = op2.Dat(nodes, np.zeros(nodes.size, dtype=valuetype), valuetype, "b") - x = op2.Dat(nodes, np.zeros(nodes.size, dtype=valuetype), valuetype, "x") - - # Set up initial condition - - f_vals = np.array([2 * X + 4 * Y for X, Y in coords.data], dtype=valuetype) - f = op2.Dat(nodes, f_vals, valuetype, "f") - - # Assemble and solve - - op2.par_loop(mass, elements, - mat(op2.INC, (elem_node[op2.i[0]], elem_node[op2.i[1]])), - coords(op2.READ, elem_node, flatten=True)) - - op2.par_loop(rhs, elements, - b(op2.INC, elem_node[op2.i[0]]), - coords(op2.READ, elem_node, flatten=True), - f(op2.READ, elem_node)) - - solver = op2.Solver() - solver.solve(mat, x, b) - - # Print solution (if necessary) - if opt['print_output']: - print "Expected solution: %s" % f.data - print "Computed solution: %s" % x.data - - # Save output (if necessary) - if opt['return_output']: - return f.data, x.data - if opt['save_output']: - from cPickle import dump, HIGHEST_PROTOCOL - import gzip - out = gzip.open("mass2d_triangle.out.gz", "wb") - dump((f.data, x.data, b.data, mat.array), out, HIGHEST_PROTOCOL) - out.close() - -parser = utils.parser(group=True, description=__doc__) -parser.add_argument('-m', '--mesh', required=True, - help='Base name of triangle mesh \ - (excluding the .ele or .node extension)') -parser.add_argument('-r', '--return-output', action='store_true', - help='Return output for testing') -parser.add_argument('-s', '--save-output', action='store_true', - help='Save the output of the run (used for testing)') -parser.add_argument('--print-output', action='store_true', - help='Print the output of the run to stdout') -parser.add_argument('-p', '--profile', action='store_true', - help='Create a cProfile for the run') - -if __name__ == '__main__': - opt = vars(parser.parse_args()) - op2.init(**opt) - - if opt['profile']: - import cProfile - filename = 'mass2d_triangle.%s.cprofile' % os.path.split(opt['mesh'])[-1] - cProfile.run('main(opt)', filename=filename) - else: - main(opt) From b3898f094a4e1ab1aa19ac545baa44313ea09be8 Mon Sep 17 00:00:00 2001 From: Florian Rathgeber Date: Fri, 4 Apr 2014 14:33:26 +0100 Subject: [PATCH 12/18] Remove vector mass demo --- demo/mass_vector_ffc.py | 137 ---------------------------------------- 1 file changed, 137 deletions(-) delete mode 100644 demo/mass_vector_ffc.py diff --git a/demo/mass_vector_ffc.py b/demo/mass_vector_ffc.py deleted file mode 100644 index e29eec9be..000000000 --- a/demo/mass_vector_ffc.py +++ /dev/null @@ -1,137 +0,0 @@ -# This file is part of PyOP2 -# -# PyOP2 is Copyright (c) 2012, Imperial College London and -# others. Please see the AUTHORS file in the main source directory for -# a full list of copyright holders. All rights reserved. -# -# Redistribution and use in source and binary forms, with or without -# modification, are permitted provided that the following conditions -# are met: -# -# * Redistributions of source code must retain the above copyright -# notice, this list of conditions and the following disclaimer. -# * Redistributions in binary form must reproduce the above copyright -# notice, this list of conditions and the following disclaimer in the -# documentation and/or other materials provided with the distribution. -# * The name of Imperial College London or that of other -# contributors may not be used to endorse or promote products -# derived from this software without specific prior written -# permission. -# -# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTERS -# ''AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT -# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS -# FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE -# COPYRIGHT HOLDERS OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, -# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES -# (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR -# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) -# HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, -# STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) -# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED -# OF THE POSSIBILITY OF SUCH DAMAGE. - -"""PyOP2 2D mass equation demo (vector field version) - -This demo solves the identity equation for a vector variable on a quadrilateral -domain. The initial condition is that all DoFs are [1, 2]^T - -This demo requires the MAPDES forks of FFC, FIAT and UFL which are found at: - - https://bitbucket.org/mapdes/ffc - https://bitbucket.org/mapdes/fiat - https://bitbucket.org/mapdes/ufl -""" - -from pyop2 import op2, utils -from ufl import * -from pyop2.ffc_interface import compile_form - -import numpy as np - - -def main(opt): - # Set up finite element identity problem - - E = VectorElement("Lagrange", "triangle", 1) - - v = TestFunction(E) - u = TrialFunction(E) - f = Coefficient(E) - - a = inner(v, u) * dx - L = inner(v, f) * dx - - # Generate code for mass and rhs assembly. - - mass, = compile_form(a, "mass") - rhs, = compile_form(L, "rhs") - - # Set up simulation data structures - - NUM_ELE = 2 - NUM_NODES = 4 - valuetype = np.float64 - - nodes = op2.Set(NUM_NODES, "nodes") - elements = op2.Set(NUM_ELE, "elements") - - elem_node_map = np.array([0, 1, 3, 2, 3, 1], dtype=np.uint32) - elem_vnode = op2.Map(elements, nodes, 3, elem_node_map, "elem_vnode") - - sparsity = op2.Sparsity(nodes ** 2, elem_vnode, "sparsity") - mat = op2.Mat(sparsity, valuetype, "mat") - - coord_vals = np.array([(0.0, 0.0), (2.0, 0.0), (1.0, 1.0), (0.0, 1.5)], - dtype=valuetype) - coords = op2.Dat(nodes ** 2, coord_vals, valuetype, "coords") - - f = op2.Dat(nodes ** 2, np.array([(1.0, 2.0)] * 4), valuetype, "f") - b = op2.Dat(nodes ** 2, np.zeros(2 * NUM_NODES), valuetype, "b") - x = op2.Dat(nodes ** 2, np.zeros(2 * NUM_NODES), valuetype, "x") - - # Assemble and solve - - op2.par_loop(mass, elements, - mat(op2.INC, (elem_vnode[op2.i[0]], elem_vnode[op2.i[1]]), flatten=True), - coords(op2.READ, elem_vnode, flatten=True)) - - op2.par_loop(rhs, elements, - b(op2.INC, elem_vnode[op2.i[0]], flatten=True), - coords(op2.READ, elem_vnode, flatten=True), - f(op2.READ, elem_vnode, flatten=True)) - - solver = op2.Solver() - solver.solve(mat, x, b) - - # Print solution - if opt['print_output']: - print "Expected solution: %s" % f.data - print "Computed solution: %s" % x.data - - # Save output (if necessary) - if opt['return_output']: - return f.data, x.data - if opt['save_output']: - import pickle - with open("mass_vector.out", "w") as out: - pickle.dump((f.data, x.data), out) - -parser = utils.parser(group=True, description=__doc__) -parser.add_argument('--print-output', action='store_true', help='Print output') -parser.add_argument('-r', '--return-output', action='store_true', - help='Return output for testing') -parser.add_argument('-s', '--save-output', action='store_true', - help='Save the output of the run (used for testing)') -parser.add_argument('-p', '--profile', action='store_true', - help='Create a cProfile for the run') - -if __name__ == '__main__': - opt = vars(parser.parse_args()) - op2.init(**opt) - - if opt['profile']: - import cProfile - cProfile.run('main(opt)', filename='mass_vector_ffc.cprofile') - else: - main(opt) From 49a2508269a8b090b6a6d59f602bec3e3a93bf61 Mon Sep 17 00:00:00 2001 From: Florian Rathgeber Date: Fri, 4 Apr 2014 14:56:24 +0100 Subject: [PATCH 13/18] Remove weak bcs demo --- demo/weak_bcs_ffc.py | 204 ------------------------------------------- 1 file changed, 204 deletions(-) delete mode 100644 demo/weak_bcs_ffc.py diff --git a/demo/weak_bcs_ffc.py b/demo/weak_bcs_ffc.py deleted file mode 100644 index 4ae4c24ed..000000000 --- a/demo/weak_bcs_ffc.py +++ /dev/null @@ -1,204 +0,0 @@ -# This file is part of PyOP2 -# -# PyOP2 is Copyright (c) 2012, Imperial College London and -# others. Please see the AUTHORS file in the main source directory for -# a full list of copyright holders. All rights reserved. -# -# Redistribution and use in source and binary forms, with or without -# modification, are permitted provided that the following conditions -# are met: -# -# * Redistributions of source code must retain the above copyright -# notice, this list of conditions and the following disclaimer. -# * Redistributions in binary form must reproduce the above copyright -# notice, this list of conditions and the following disclaimer in the -# documentation and/or other materials provided with the distribution. -# * The name of Imperial College London or that of other -# contributors may not be used to endorse or promote products -# derived from this software without specific prior written -# permission. -# -# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTERS -# ''AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT -# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS -# FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE -# COPYRIGHT HOLDERS OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, -# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES -# (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR -# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) -# HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, -# STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) -# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED -# OF THE POSSIBILITY OF SUCH DAMAGE. - -"""PyOP2 laplace equation demo (weak BCs) - -This demo uses ffc-generated kernels to solve the Laplace equation on a unit -square with boundary conditions: - - u = 1 on y = 0 - du/dn = 2 on y = 1 - -The domain is meshed as follows: - - *-*-* - |/|/| - *-*-* - |/|/| - *-*-* - -This demo requires the MAPDES forks of FFC, FIAT and UFL which are found at: - - https://bitbucket.org/mapdes/ffc - https://bitbucket.org/mapdes/fiat - https://bitbucket.org/mapdes/ufl -""" - -from pyop2 import op2, utils -from pyop2.ffc_interface import compile_form -from ufl import * - -import numpy as np - - -def main(opt): - # Set up finite element problem - - E = FiniteElement("Lagrange", "triangle", 1) - - v = TestFunction(E) - u = TrialFunction(E) - f = Coefficient(E) - g = Coefficient(E) - - a = dot(grad(v,), grad(u)) * dx - L = v * f * dx + v * g * ds(2) - - # Generate code for Laplacian and rhs assembly. - - laplacian, = compile_form(a, "laplacian") - rhs, weak = compile_form(L, "rhs") - - # Set up simulation data structures - - NUM_ELE = 8 - NUM_NODES = 9 - NUM_BDRY_ELE = 2 - NUM_BDRY_NODE = 3 - valuetype = np.float64 - - nodes = op2.Set(NUM_NODES, "nodes") - elements = op2.Set(NUM_ELE, "elements") - - # Elements that Weak BC will be assembled over - top_bdry_elements = op2.Set(NUM_BDRY_ELE, "top_boundary_elements") - # Nodes that Strong BC will be applied over - bdry_nodes = op2.Set(NUM_BDRY_NODE, "boundary_nodes") - - elem_node_map = np.array([0, 1, 4, 4, 3, 0, 1, 2, 5, 5, 4, 1, 3, 4, 7, 7, - 6, 3, 4, 5, 8, 8, 7, 4], dtype=np.uint32) - elem_node = op2.Map(elements, nodes, 3, elem_node_map, "elem_node") - - top_bdry_elem_node_map = np.array([7, 6, 3, 8, 7, 4], dtype=valuetype) - top_bdry_elem_node = op2.Map(top_bdry_elements, nodes, 3, - top_bdry_elem_node_map, "top_bdry_elem_node") - - bdry_node_node_map = np.array([0, 1, 2], dtype=valuetype) - bdry_node_node = op2.Map( - bdry_nodes, nodes, 1, bdry_node_node_map, "bdry_node_node") - - sparsity = op2.Sparsity((nodes, nodes), (elem_node, elem_node), "sparsity") - mat = op2.Mat(sparsity, valuetype, "mat") - - coord_vals = np.array([(0.0, 0.0), (0.5, 0.0), (1.0, 0.0), - (0.0, 0.5), (0.5, 0.5), (1.0, 0.5), - (0.0, 1.0), (0.5, 1.0), (1.0, 1.0)], - dtype=valuetype) - coords = op2.Dat(nodes ** 2, coord_vals, valuetype, "coords") - - u_vals = np.array([1.0, 1.0, 1.0, 2.0, 2.0, 2.0, 3.0, 3.0, 3.0]) - f = op2.Dat(nodes, np.zeros(NUM_NODES, dtype=valuetype), valuetype, "f") - b = op2.Dat(nodes, np.zeros(NUM_NODES, dtype=valuetype), valuetype, "b") - x = op2.Dat(nodes, np.zeros(NUM_NODES, dtype=valuetype), valuetype, "x") - u = op2.Dat(nodes, u_vals, valuetype, "u") - - bdry = op2.Dat(bdry_nodes, np.ones(3, dtype=valuetype), valuetype, "bdry") - - # This isn't perfect, defining the boundary gradient on more nodes than are on - # the boundary is couter-intuitive - bdry_grad_vals = np.asarray([2.0] * 9, dtype=valuetype) - bdry_grad = op2.Dat(nodes, bdry_grad_vals, valuetype, "gradient") - facet = op2.Global(1, 2, np.uint32, "facet") - - # If a form contains multiple integrals with differing coefficients, FFC - # generates kernels that take all the coefficients of the entire form (not - # only the respective integral) as arguments. Arguments that correspond to - # forms that are not used in that integral are simply not referenced. - # We therefore need a dummy argument in place of the coefficient that is not - # used in the par_loop for OP2 to generate the correct kernel call. - - # Assemble matrix and rhs - - op2.par_loop(laplacian, elements, - mat(op2.INC, (elem_node[op2.i[0]], elem_node[op2.i[1]])), - coords(op2.READ, elem_node, flatten=True)) - - op2.par_loop(rhs, elements, - b(op2.INC, elem_node[op2.i[0]]), - coords(op2.READ, elem_node, flatten=True), - f(op2.READ, elem_node), - bdry_grad(op2.READ, elem_node)) # argument ignored - - # Apply weak BC - - op2.par_loop(weak, top_bdry_elements, - b(op2.INC, top_bdry_elem_node[op2.i[0]]), - coords(op2.READ, top_bdry_elem_node, flatten=True), - f(op2.READ, top_bdry_elem_node), # argument ignored - bdry_grad(op2.READ, top_bdry_elem_node), - facet(op2.READ)) - - # Apply strong BC - - mat.zero_rows([0, 1, 2], 1.0) - strongbc_rhs = op2.Kernel(""" - void strongbc_rhs(double *val, double *target) { *target = *val; } - """, "strongbc_rhs") - op2.par_loop(strongbc_rhs, bdry_nodes, - bdry(op2.READ), - b(op2.WRITE, bdry_node_node[0])) - - solver = op2.Solver(ksp_type='gmres') - solver.solve(mat, x, b) - - # Print solution - if opt['return_output']: - return u.data, x.data - if opt['print_output']: - print "Expected solution: %s" % u.data - print "Computed solution: %s" % x.data - - # Save output (if necessary) - if opt['save_output']: - import pickle - with open("weak_bcs.out", "w") as out: - pickle.dump((u.data, x.data), out) - -parser = utils.parser(group=True, description=__doc__) -parser.add_argument('--print-output', action='store_true', help='Print output') -parser.add_argument('-r', '--return-output', action='store_true', - help='Return output for testing') -parser.add_argument('-s', '--save-output', action='store_true', - help='Save the output of the run (used for testing)') -parser.add_argument('-p', '--profile', action='store_true', - help='Create a cProfile for the run') - -if __name__ == '__main__': - opt = vars(parser.parse_args()) - op2.init(**opt) - - if opt['profile']: - import cProfile - cProfile.run('main(opt)', filename='weak_bcs_ffc.cprofile') - else: - main(opt) From e0c1c57d0d8b91dbd16a3ac84ec56ea531fabe8c Mon Sep 17 00:00:00 2001 From: Florian Rathgeber Date: Thu, 30 Jan 2014 17:57:49 +0000 Subject: [PATCH 14/18] Remove FFC interface --- pyop2/__init__.py | 1 - pyop2/ffc_interface.py | 143 ----------------- pyop2/pyop2_geometry.h | 274 -------------------------------- setup.py | 2 +- test/unit/test_ffc_interface.py | 133 ---------------- 5 files changed, 1 insertion(+), 552 deletions(-) delete mode 100644 pyop2/ffc_interface.py delete mode 100644 pyop2/pyop2_geometry.h delete mode 100644 test/unit/test_ffc_interface.py diff --git a/pyop2/__init__.py b/pyop2/__init__.py index a5e8a5088..e4ba94aaf 100644 --- a/pyop2/__init__.py +++ b/pyop2/__init__.py @@ -8,4 +8,3 @@ from op2 import * from version import __version__, __version_info__ # noqa: we just want to expose these -from ffc_interface import clear_cache as clear_ffc_cache # noqa: expose to user diff --git a/pyop2/ffc_interface.py b/pyop2/ffc_interface.py deleted file mode 100644 index b13c7ae6d..000000000 --- a/pyop2/ffc_interface.py +++ /dev/null @@ -1,143 +0,0 @@ -# This file is part of PyOP2 -# -# PyOP2 is Copyright (c) 2012, Imperial College London and -# others. Please see the AUTHORS file in the main source directory for -# a full list of copyright holders. All rights reserved. -# -# Redistribution and use in source and binary forms, with or without -# modification, are permitted provided that the following conditions -# are met: -# -# * Redistributions of source code must retain the above copyright -# notice, this list of conditions and the following disclaimer. -# * Redistributions in binary form must reproduce the above copyright -# notice, this list of conditions and the following disclaimer in the -# documentation and/or other materials provided with the distribution. -# * The name of Imperial College London or that of other -# contributors may not be used to endorse or promote products -# derived from this software without specific prior written -# permission. -# -# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTERS -# ''AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT -# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS -# FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE -# COPYRIGHT HOLDERS OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, -# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES -# (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR -# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) -# HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, -# STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) -# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED -# OF THE POSSIBILITY OF SUCH DAMAGE. - -"""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 ffc.log import set_level, ERROR - -from caching import DiskCached -from op2 import Kernel -from mpi import MPI - -from ir.ast_base import PreprocessNode, Root - -_form_cache = {} - -# Silence FFC -set_level(ERROR) - -ffc_parameters = default_parameters() -ffc_parameters['write_file'] = False -ffc_parameters['format'] = 'pyop2' -ffc_parameters['pyop2-ir'] = True - -# Include an md5 hash of pyop2_geometry.h in the cache key -with open(os.path.join(os.path.dirname(__file__), 'pyop2_geometry.h')) as f: - _pyop2_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'))) - - -class FFCKernel(DiskCached): - - _cache = {} - _cachedir = os.path.join(tempfile.gettempdir(), - 'pyop2-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__ + - _pyop2_geometry_md5 + constants.FFC_VERSION + - constants.PYOP2_VERSION).hexdigest() - - def __init__(self, form, name): - if self._initialized: - return - - incl = PreprocessNode('#include "pyop2_geometry.h"\n') - 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)) - self.kernels = tuple(kernels) - - self._initialized = True - - -def compile_form(form, name): - """Compile a form using FFC and return a :class:`pyop2.op2.Kernel`.""" - - # 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/pyop2/pyop2_geometry.h b/pyop2/pyop2_geometry.h deleted file mode 100644 index 5ef324927..000000000 --- a/pyop2/pyop2_geometry.h +++ /dev/null @@ -1,274 +0,0 @@ -/* --- 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/setup.py b/setup.py index a1d9c6f00..0d8473a71 100644 --- a/setup.py +++ b/setup.py @@ -134,7 +134,7 @@ def run(self): test_requires=test_requires, packages=['pyop2', 'pyop2.ir', 'pyop2_utils'], package_data={ - 'pyop2': ['assets/*', 'mat_utils.*', 'sparsity_utils.*', '*.pyx', 'pyop2_geometry.h']}, + 'pyop2': ['assets/*', 'mat_utils.*', 'sparsity_utils.*', '*.pyx']}, scripts=glob('scripts/*'), cmdclass=cmdclass, ext_modules=[NumpyExtension('pyop2.plan', plan_sources), diff --git a/test/unit/test_ffc_interface.py b/test/unit/test_ffc_interface.py deleted file mode 100644 index 3119b7afb..000000000 --- a/test/unit/test_ffc_interface.py +++ /dev/null @@ -1,133 +0,0 @@ -# This file is part of PyOP2 -# -# PyOP2 is Copyright (c) 2012, Imperial College London and -# others. Please see the AUTHORS file in the main source directory for -# a full list of copyright holders. All rights reserved. -# -# Redistribution and use in source and binary forms, with or without -# modification, are permitted provided that the following conditions -# are met: -# -# * Redistributions of source code must retain the above copyright -# notice, this list of conditions and the following disclaimer. -# * Redistributions in binary form must reproduce the above copyright -# notice, this list of conditions and the following disclaimer in the -# documentation and/or other materials provided with the distribution. -# * The name of Imperial College London or that of other -# contributors may not be used to endorse or promote products -# derived from this software without specific prior written -# permission. -# -# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTERS -# ''AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT -# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS -# FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE -# COPYRIGHT HOLDERS OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, -# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES -# (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR -# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) -# HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, -# STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) -# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED -# OF THE POSSIBILITY OF SUCH DAMAGE. - -import pytest -ffc_interface = pytest.importorskip('pyop2.ffc_interface') -import os -from ufl import * - - -@pytest.mark.xfail("not hasattr(ffc_interface.constants, 'PYOP2_VERSION')") -class TestFFCCache: - - """FFC code generation cache tests.""" - - @pytest.fixture - def mass(cls): - e = FiniteElement('CG', triangle, 1) - u = TestFunction(e) - v = TrialFunction(e) - return u * v * dx - - @pytest.fixture - def mass2(cls): - e = FiniteElement('CG', triangle, 2) - u = TestFunction(e) - v = TrialFunction(e) - return u * v * dx - - @pytest.fixture - def rhs(cls): - e = FiniteElement('CG', triangle, 1) - v = TrialFunction(e) - g = Coefficient(e) - return g * v * ds - - @pytest.fixture - def rhs2(cls): - e = FiniteElement('CG', triangle, 1) - v = TrialFunction(e) - f = Coefficient(e) - g = Coefficient(e) - return f * v * dx + g * v * ds - - @pytest.fixture - def cache_key(cls, mass): - return ffc_interface.FFCKernel(mass, 'mass').cache_key - - def test_ffc_cache_dir_exists(self, backend): - """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, backend, 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, backend, 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, backend, mass): - """Compiling a form attaches form data.""" - ffc_interface.compile_form(mass, 'mass') - - assert mass.form_data() - - def test_ffc_same_form(self, backend, 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, backend, mass, mass2): - """Compiling different forms should not load kernels from cache.""" - k1 = ffc_interface.compile_form(mass, 'mass') - k2 = ffc_interface.compile_form(mass2, 'mass') - - assert k1 is not k2 - - def test_ffc_different_names(self, backend, mass): - """Compiling different forms should not load kernels from cache.""" - k1 = ffc_interface.compile_form(mass, 'mass') - k2 = ffc_interface.compile_form(mass, 'mass2') - - assert k1 is not k2 - - def test_ffc_cell_kernel(self, backend, 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, backend, 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, backend, 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 eb4ce0e44b2861160cb17122b57286a3203e1af5 Mon Sep 17 00:00:00 2001 From: Florian Rathgeber Date: Fri, 31 Jan 2014 10:18:55 +0000 Subject: [PATCH 15/18] Remove compatible FFC version --- pyop2/version.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/pyop2/version.py b/pyop2/version.py index 562e055cc..2bfaa0eba 100644 --- a/pyop2/version.py +++ b/pyop2/version.py @@ -1,4 +1,2 @@ __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__)) From f0935e1f682892ca63af13a0846d823f7ae89e2d Mon Sep 17 00:00:00 2001 From: Lawrence Mitchell Date: Tue, 25 Mar 2014 18:43:19 +0000 Subject: [PATCH 16/18] Repurpose pyop2-clean to remove cached compiled libraries Now that the ffc interface is gone, pyop2-clean doesn't need to remove cached kernels anymore. However, we do need a way to blow away compiled libraries which we now cache ourselves. --- pyop2/compilation.py | 34 ++++++++++++++++++++++++++++++++++ scripts/pyop2-clean | 6 ++---- 2 files changed, 36 insertions(+), 4 deletions(-) diff --git a/pyop2/compilation.py b/pyop2/compilation.py index 000aee212..0d27c2357 100644 --- a/pyop2/compilation.py +++ b/pyop2/compilation.py @@ -203,3 +203,37 @@ def load(src, fn_name, cppargs=[], ldargs=[], argtypes=None, restype=None): fn.argtypes = argtypes fn.restype = restype return fn + + +def clear_cache(prompt=False): + """Clear the PyOP2 compiler cache. + + :arg prompt: if ``True`` prompt before removing any files + """ + cachedir = configuration['cache_dir'] + + files = [os.path.join(cachedir, f) for f in os.listdir(cachedir) + if os.path.isfile(os.path.join(cachedir, f))] + nfiles = len(files) + + if nfiles == 0: + print "No cached libraries to remove" + return + + remove = True + if prompt: + + user = raw_input("Remove %d cached libraries from %s? [Y/n]: " % (nfiles, cachedir)) + + while user.lower() not in ['', 'y', 'n']: + print "Please answer y or n." + user = raw_input("Remove %d cached libraries from %s? [Y/n]: " % (nfiles, cachedir)) + + if user.lower() == 'n': + remove = False + + if remove: + print "Removing %d cached libraries from %s" % (nfiles, cachedir) + [os.remove(f) for f in files] + else: + print "Not removing cached libraries" diff --git a/scripts/pyop2-clean b/scripts/pyop2-clean index 931e0e5cc..ab29f1245 100755 --- a/scripts/pyop2-clean +++ b/scripts/pyop2-clean @@ -1,8 +1,6 @@ #!/usr/bin/env python - -from pyop2.ffc_interface import clear_cache, FFCKernel +from pyop2.compilation import clear_cache if __name__ == '__main__': - print 'Removing cached ffc kernels from %s' % FFCKernel._cachedir - clear_cache() + clear_cache(prompt=True) From ec26be209a3e85a49af4f3c52620f3fab2680d0d Mon Sep 17 00:00:00 2001 From: Florian Rathgeber Date: Fri, 4 Apr 2014 15:23:12 +0100 Subject: [PATCH 17/18] README: remove FFC, FIAT, UFL dependencies --- README.rst | 24 ------------------------ 1 file changed, 24 deletions(-) diff --git a/README.rst b/README.rst index f2e2b33bb..37bcf8616 100644 --- a/README.rst +++ b/README.rst @@ -320,21 +320,6 @@ When installing PyOP2 via ``python setup.py install`` the extension modules will be built automatically and amending ``$PYTHONPATH`` is not necessary. -FFC Interface -------------- - -Solving UFL_ finite element equations requires a fork of FFC_, UFL_ -and FIAT_. Note that FFC_ requires a version of Instant_. - -Install FFC_ and all dependencies via pip:: - - sudo pip install \ - git+https://bitbucket.org/mapdes/ffc.git#egg=ffc - git+https://bitbucket.org/mapdes/ufl.git#egg=ufl - git+https://bitbucket.org/mapdes/fiat.git#egg=fiat - git+https://bitbucket.org/fenics-project/instant.git#egg=instant - hg+https://bitbucket.org/khinsen/scientificpython - Setting up the environment -------------------------- @@ -347,12 +332,6 @@ definitions as necessary:: export PETSC_DIR=/path/to/petsc export PETSC_ARCH=linux-gnu-c-opt - #Add UFL and FFC to PYTHONPATH if in non-standard location - export UFL_DIR=/path/to/ufl - export FFC_DIR=/path/to/ffc - export PYTHONPATH=$UFL_DIR:$FFC_DIR:$PYTHONPATH - # Add any other Python module in non-standard locations - #Add PyOP2 to PYTHONPATH export PYTHONPATH=/path/to/PyOP2:$PYTHONPATH @@ -417,7 +396,4 @@ manner as required. .. _PPA: https://launchpad.net/~amcg/+archive/petsc3.4/ .. _PETSc: http://www.mcs.anl.gov/petsc/ .. _petsc4py: http://pythonhosted.org/petsc4py/ -.. _FFC: https://bitbucket.org/mapdes/ffc -.. _FIAT: https://bitbucket.org/mapdes/fiat -.. _UFL: https://bitbucket.org/mapdes/ufl .. _Instant: https://bitbucket.org/fenics-project/instant From 52bc27d676623f43b5ee6713dd78f8cc1fcc4830 Mon Sep 17 00:00:00 2001 From: Florian Rathgeber Date: Fri, 4 Apr 2014 15:24:52 +0100 Subject: [PATCH 18/18] Remove FFC dependencies from requirements and install.sh --- install.sh | 10 ---------- requirements-minimal.txt | 5 ----- 2 files changed, 15 deletions(-) diff --git a/install.sh b/install.sh index 2c62d6427..618a3b8cf 100644 --- a/install.sh +++ b/install.sh @@ -64,16 +64,6 @@ echo | tee -a $LOGFILE # Install Cython so we can build PyOP2 from source ${PIP} Cython decorator numpy >> $LOGFILE 2>&1 -echo "*** Installing FEniCS dependencies ***" | tee -a $LOGFILE -echo | tee -a $LOGFILE - -${PIP} \ - git+https://bitbucket.org/mapdes/ffc#egg=ffc \ - git+https://bitbucket.org/mapdes/ufl#egg=ufl \ - git+https://bitbucket.org/mapdes/fiat#egg=fiat \ - git+https://bitbucket.org/fenics-project/instant#egg=instant \ - hg+https://bitbucket.org/khinsen/scientificpython >> $LOGFILE 2>&1 - echo "*** Installing PETSc ***" | tee -a $LOGFILE echo | tee -a $LOGFILE diff --git a/requirements-minimal.txt b/requirements-minimal.txt index 29a33b4c3..66039956a 100644 --- a/requirements-minimal.txt +++ b/requirements-minimal.txt @@ -11,11 +11,6 @@ pytest>=2.3 flake8>=2.1.0 pycparser>=2.10 mpi4py>=1.3.1 -git+https://bitbucket.org/fenics-project/instant.git#egg=instant -git+https://bitbucket.org/mapdes/ufl.git#egg=ufl -git+https://bitbucket.org/mapdes/fiat.git#egg=fiat -git+https://bitbucket.org/mapdes/ffc.git#egg=ffc -hg+https://bitbucket.org/khinsen/scientificpython h5py>=2.0.0 git+https://bitbucket.org/petsc/petsc.git#egg=petsc git+https://bitbucket.org/petsc/petsc4py.git#egg=petsc4py