tsfc/tests/test_underintegration.py at 5935641d4cf6c9f50e1e11d3276097ea37841ba8 · firedrakeproject/tsfc · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
from functools import reduce

import numpy
import pytest

from coffee.visitors import EstimateFlops

from ufl import (Mesh, FunctionSpace, FiniteElement, VectorElement,
                 TestFunction, TrialFunction, TensorProductCell, dx,
                 action, interval, quadrilateral, dot, grad)

from FIAT import ufc_cell
from FIAT.quadrature import GaussLobattoLegendreQuadratureLineRule, GaussLegendreQuadratureLineRule, ExtendedGaussLegendreQuadratureLineRule

from finat.point_set import GaussLobattoLegendrePointSet, GaussLegendrePointSet, ExtendedGaussLegendrePointSet
from finat.quadrature import QuadratureRule, TensorProductQuadratureRule

from tsfc import compile_form


def gll_quadrature_rule(cell, elem_deg):
    fiat_cell = ufc_cell("interval")
    fiat_rule = GaussLobattoLegendreQuadratureLineRule(fiat_cell, elem_deg + 1)
    line_rules = [QuadratureRule(GaussLobattoLegendrePointSet(fiat_rule.get_points()),
                                 fiat_rule.get_weights())
                  for _ in range(cell.topological_dimension())]
    finat_rule = reduce(lambda a, b: TensorProductQuadratureRule([a, b]), line_rules)
    return finat_rule


def gl_quadrature_rule(cell, elem_deg):
    fiat_cell = ufc_cell("interval")
    fiat_rule = GaussLegendreQuadratureLineRule(fiat_cell, elem_deg + 1)
    line_rules = [QuadratureRule(GaussLegendrePointSet(fiat_rule.get_points()),
                                 fiat_rule.get_weights())
                  for _ in range(cell.topological_dimension())]
    finat_rule = reduce(lambda a, b: TensorProductQuadratureRule([a, b]), line_rules)
    return finat_rule


def egl_quadrature_rule(cell, elem_deg):
    fiat_cell = ufc_cell("interval")
    fiat_rule = ExtendedGaussLegendreQuadratureLineRule(fiat_cell, elem_deg + 1)
    line_rules = [QuadratureRule(ExtendedGaussLegendrePointSet(fiat_rule.get_points()),
                                 fiat_rule.get_weights())
                  for _ in range(cell.topological_dimension())]
    finat_rule = reduce(lambda a, b: TensorProductQuadratureRule([a, b]), line_rules)
    return finat_rule


def mass_cg(cell, degree):
    m = Mesh(VectorElement('Q', cell, 1))
    V = FunctionSpace(m, FiniteElement('Q', cell, degree, variant='spectral'))
    u = TrialFunction(V)
    v = TestFunction(V)
    return u*v*dx(scheme=gll_quadrature_rule(cell, degree))


def mass_cg_egl(cell, degree):
    m = Mesh(VectorElement('Q', cell, 1))
    V = FunctionSpace(m, FiniteElement('Q', cell, degree, variant='dualmse'))
    u = TrialFunction(V)
    v = TestFunction(V)
    return u*v*dx(scheme=egl_quadrature_rule(cell, degree))


def mass_dg(cell, degree):
    m = Mesh(VectorElement('Q', cell, 1))
    V = FunctionSpace(m, FiniteElement('DQ', cell, degree, variant='spectral'))
    u = TrialFunction(V)
    v = TestFunction(V)
    return u*v*dx(scheme=gl_quadrature_rule(cell, degree))


def laplace(cell, degree):
    m = Mesh(VectorElement('Q', cell, 1))
    V = FunctionSpace(m, FiniteElement('Q', cell, degree, variant='spectral'))
    u = TrialFunction(V)
    v = TestFunction(V)
    return dot(grad(u), grad(v))*dx(scheme=gll_quadrature_rule(cell, degree))


def laplace_dualmse(cell, degree):
    m = Mesh(VectorElement('Q', cell, 1))
    V = FunctionSpace(m, FiniteElement('Q', cell, degree, variant='dualmse'))
    u = TrialFunction(V)
    v = TestFunction(V)
    return dot(grad(u), grad(v))*dx(scheme=egl_quadrature_rule(cell, degree))


def count_flops(form):
    kernel, = compile_form(form, parameters=dict(mode='spectral'))
    return EstimateFlops().visit(kernel.ast)


@pytest.mark.parametrize('form', [mass_cg, mass_dg, mass_cg_egl])
@pytest.mark.parametrize(('cell', 'order'),
                         [(quadrilateral, 2),
                          (TensorProductCell(interval, interval), 2),
                          (TensorProductCell(quadrilateral, interval), 3)])
def test_mass(form, cell, order):
    degrees = numpy.arange(4, 10)
    flops = [count_flops(form(cell, int(degree))) for degree in degrees]
    rates = numpy.diff(numpy.log(flops)) / numpy.diff(numpy.log(degrees + 1))
    assert (rates < order).all()


@pytest.mark.parametrize('form', [mass_cg, mass_dg, mass_cg_egl])
@pytest.mark.parametrize(('cell', 'order'),
                         [(quadrilateral, 2),
                          (TensorProductCell(interval, interval), 2),
                          (TensorProductCell(quadrilateral, interval), 3)])
def test_mass_action(form, cell, order):
    degrees = numpy.arange(4, 10)
    flops = [count_flops(action(form(cell, int(degree)))) for degree in degrees]
    rates = numpy.diff(numpy.log(flops)) / numpy.diff(numpy.log(degrees + 1))
    assert (rates < order).all()


@pytest.mark.parametrize(('cell', 'order'),
                         [(quadrilateral, 4),
                          (TensorProductCell(interval, interval), 4),
                          (TensorProductCell(quadrilateral, interval), 5)])
def test_laplace(cell, order):
    degrees = numpy.arange(4, 10)
    flops = [count_flops(laplace(cell, int(degree))) for degree in degrees]
    rates = numpy.diff(numpy.log(flops)) / numpy.diff(numpy.log(degrees + 1))
    assert (rates < order).all()


@pytest.mark.parametrize(('cell', 'order'),
                         [(quadrilateral, 4),
                          (TensorProductCell(interval, interval), 4),
                          (TensorProductCell(quadrilateral, interval), 5)])
def test_laplace_dualmse(cell, order):
    degrees = numpy.arange(4, 10)
    flops = [count_flops(laplace_dualmse(cell, int(degree))) for degree in degrees]
    rates = numpy.diff(numpy.log(flops)) / numpy.diff(numpy.log(degrees + 1))
    assert (rates < order).all()


if __name__ == "__main__":
    import os
    import sys
    pytest.main(args=[os.path.abspath(__file__)] + sys.argv[1:])