Test Alfeld-Sorokina #5927
8077 tests run, 7279 passed, 786 skipped, 12 failed.
Annotations
Check failure on line 57 in tests/regression/test_facet_split.py
github-actions / Firedrake real
test_facet_split.test_facet_split[lu-True]
NotImplementedError: unsupported functional type
Raw output
quadrilateral = True, pc_type = 'lu'
@pytest.mark.parametrize("quadrilateral", [True, False])
@pytest.mark.parametrize("pc_type", ["lu", "jacobi"])
def test_facet_split(quadrilateral, pc_type):
> assert run_facet_split(quadrilateral, pc_type) < 1E-10
tests/regression/test_facet_split.py:64:
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
tests/regression/test_facet_split.py:57: in run_facet_split
solve(a == L, uh, bcs=bcs, solver_parameters=parameters)
petsc4py/PETSc/Log.pyx:188: in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
???
petsc4py/PETSc/Log.pyx:189: in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
???
firedrake/adjoint_utils/solving.py:57: in wrapper
output = solve(*args, **kwargs)
firedrake/solving.py:140: in solve
_solve_varproblem(*args, **kwargs)
firedrake/solving.py:177: in _solve_varproblem
solver.solve()
petsc4py/PETSc/Log.pyx:188: in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
???
petsc4py/PETSc/Log.pyx:189: in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
???
firedrake/adjoint_utils/variational_solver.py:89: in wrapper
out = solve(self, **kwargs)
firedrake/variational_solver.py:307: in solve
dbc.apply(problem.u_restrict)
petsc4py/PETSc/Log.pyx:188: in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
???
petsc4py/PETSc/Log.pyx:189: in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
???
firedrake/adjoint_utils/dirichletbc.py:32: in wrapper
ret = apply(self, *args, **kwargs)
firedrake/bcs.py:451: in apply
r.assign(self.function_arg, subset=self.node_set)
../firedrake_venv/src/PyOP2/pyop2/utils.py:62: in __get__
obj.__dict__[self.__name__] = result = self.fget(obj)
firedrake/bcs.py:179: in node_set
return op2.Subset(self._function_space.node_set, self.nodes)
../firedrake_venv/src/PyOP2/pyop2/utils.py:62: in __get__
obj.__dict__[self.__name__] = result = self.fget(obj)
firedrake/bcs.py:169: in nodes
bcnodes1.append(hermite_stride(self._function_space.boundary_nodes(ss)))
firedrake/bcs.py:143: in hermite_stride
deriv_nodes = [i for i, node in enumerate(fe.fiat_equivalent.dual) if len(node.deriv_dict) != 0]
../firedrake_venv/src/tsfc/gem/utils.py:19: in __get__
obj.__dict__[self.__name__] = result = self.fget(obj)
../firedrake_venv/src/FInAT/finat/cube.py:65: in fiat_equivalent
return FIAT_FlattenedDimensions(self.product.fiat_equivalent)
../firedrake_venv/src/tsfc/gem/utils.py:19: in __get__
obj.__dict__[self.__name__] = result = self.fget(obj)
../firedrake_venv/src/FInAT/finat/tensor_product.py:82: in fiat_equivalent
return FIAT.TensorProductElement(A.fiat_equivalent, B.fiat_equivalent)
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
self = <FIAT.tensor_product.TensorProductElement object at 0x7fdfb4bb6ab0>
A = <FIAT.fdm_element.FDMLagrange object at 0x7fdfb3469160>
B = <FIAT.fdm_element.FDMLagrange object at 0x7fdfb3469160>
def __init__(self, A, B):
# set up simple things
order = min(A.get_order(), B.get_order())
if A.get_formdegree() is None or B.get_formdegree() is None:
formdegree = None
else:
formdegree = A.get_formdegree() + B.get_formdegree()
# set up reference element
ref_el = TensorProductCell(A.get_reference_element(),
B.get_reference_element())
if A.mapping()[0] != "affine" and B.mapping()[0] == "affine":
mapping = A.mapping()[0]
elif B.mapping()[0] != "affine" and A.mapping()[0] == "affine":
mapping = B.mapping()[0]
elif A.mapping()[0] == "affine" and B.mapping()[0] == "affine":
mapping = "affine"
else:
raise ValueError("check tensor product mappings - at least one must be affine")
# set up entity_ids
Adofs = A.entity_dofs()
Bdofs = B.entity_dofs()
Bsdim = B.space_dimension()
entity_ids = {}
for curAdim in Adofs:
for curBdim in Bdofs:
entity_ids[(curAdim, curBdim)] = {}
dim_cur = 0
for entityA in Adofs[curAdim]:
for entityB in Bdofs[curBdim]:
entity_ids[(curAdim, curBdim)][dim_cur] = \
[x*Bsdim + y for x in Adofs[curAdim][entityA]
for y in Bdofs[curBdim][entityB]]
dim_cur += 1
# set up dual basis
Anodes = A.dual_basis()
Bnodes = B.dual_basis()
# build the dual set by inspecting the current dual
# sets item by item.
# Currently supported cases:
# PointEval x PointEval = PointEval [scalar x scalar = scalar]
# PointScaledNormalEval x PointEval = PointScaledNormalEval [vector x scalar = vector]
# ComponentPointEvaluation x PointEval [vector x scalar = vector]
nodes = []
for Anode in Anodes:
if isinstance(Anode, functional.PointEvaluation):
for Bnode in Bnodes:
if isinstance(Bnode, functional.PointEvaluation):
# case: PointEval x PointEval
# the PointEval functional just requires the
# coordinates. these are currently stored as
# the key of a one-item dictionary. we retrieve
# these by calling get_point_dict(), and
# use the concatenation to make a new PointEval
nodes.append(functional.PointEvaluation(ref_el, _first_point(Anode) + _first_point(Bnode)))
elif isinstance(Bnode, functional.IntegralMoment):
# dummy functional for product with integral moments
nodes.append(functional.Functional(None, None, None,
{}, "Undefined"))
elif isinstance(Bnode, functional.PointDerivative):
# dummy functional for product with point derivative
nodes.append(functional.Functional(None, None, None,
{}, "Undefined"))
else:
raise NotImplementedError("unsupported functional type")
elif isinstance(Anode, functional.PointScaledNormalEvaluation):
for Bnode in Bnodes:
if isinstance(Bnode, functional.PointEvaluation):
# case: PointScaledNormalEval x PointEval
# this could be wrong if the second shape
# has spatial dimension >1, since we are not
# explicitly scaling by facet size
if len(_first_point(Bnode)) > 1:
# TODO: support this case one day
raise NotImplementedError("PointScaledNormalEval x PointEval is not yet supported if the second shape has dimension > 1")
# We cannot make a new functional.PSNEval in
# the natural way, since it tries to compute
# the normal vector by itself.
# Instead, we create things manually, and
# call Functional() with these arguments
sd = ref_el.get_spatial_dimension()
# The pt_dict is a one-item dictionary containing
# the details of the functional.
# The key is the spatial coordinate, which
# is just a concatenation of the two parts.
# The value is a list of tuples, representing
# the normal vector (scaled by the volume of
# the facet) at that point.
# Each tuple looks like (foo, (i,)); the i'th
# component of the scaled normal is foo.
# The following line is only valid when the second
# shape has spatial dimension 1 (enforced above)
Apoint, Avalue = _first_point_pair(Anode)
pt_dict = {Apoint + _first_point(Bnode): Avalue + [(0.0, (len(Apoint),))]}
# The following line should be used in the
# general case
# pt_dict = {Anode.get_point_dict().keys()[0] + Bnode.get_point_dict().keys()[0]: Anode.get_point_dict().values()[0] + [(0.0, (ii,)) for ii in range(len(Anode.get_point_dict().keys()[0]), len(Anode.get_point_dict().keys()[0]) + len(Bnode.get_point_dict().keys()[0]))]}
# THE FOLLOWING IS PROBABLY CORRECT BUT UNTESTED
shp = (sd,)
nodes.append(functional.Functional(ref_el, shp, pt_dict, {}, "PointScaledNormalEval"))
else:
raise NotImplementedError("unsupported functional type")
elif isinstance(Anode, functional.PointEdgeTangentEvaluation):
for Bnode in Bnodes:
if isinstance(Bnode, functional.PointEvaluation):
# case: PointEdgeTangentEval x PointEval
# this is very similar to the case above, so comments omitted
if len(_first_point(Bnode)) > 1:
raise NotImplementedError("PointEdgeTangentEval x PointEval is not yet supported if the second shape has dimension > 1")
sd = ref_el.get_spatial_dimension()
Apoint, Avalue = _first_point_pair(Anode)
pt_dict = {Apoint + _first_point(Bnode): Avalue + [(0.0, (len(Apoint),))]}
# THE FOLLOWING IS PROBABLY CORRECT BUT UNTESTED
shp = (sd,)
nodes.append(functional.Functional(ref_el, shp, pt_dict, {}, "PointEdgeTangent"))
else:
raise NotImplementedError("unsupported functional type")
elif isinstance(Anode, functional.ComponentPointEvaluation):
for Bnode in Bnodes:
if isinstance(Bnode, functional.PointEvaluation):
# case: ComponentPointEval x PointEval
# the CptPointEval functional requires the component
# and the coordinates. very similar to PE x PE case.
sd = ref_el.get_spatial_dimension()
nodes.append(functional.ComponentPointEvaluation(ref_el, Anode.comp, (sd,), _first_point(Anode) + _first_point(Bnode)))
else:
raise NotImplementedError("unsupported functional type")
elif isinstance(Anode, functional.FrobeniusIntegralMoment):
for Bnode in Bnodes:
if isinstance(Bnode, functional.PointEvaluation):
# case: FroIntMom x PointEval
sd = ref_el.get_spatial_dimension()
pt_dict = {}
pt_old = Anode.get_point_dict()
for pt in pt_old:
pt_dict[pt+_first_point(Bnode)] = pt_old[pt] + [(0.0, sd-1)]
# THE FOLLOWING IS PROBABLY CORRECT BUT UNTESTED
shp = (sd,)
nodes.append(functional.Functional(ref_el, shp, pt_dict, {}, "FrobeniusIntegralMoment"))
else:
raise NotImplementedError("unsupported functional type")
elif isinstance(Anode, functional.IntegralMoment):
for Bnode in Bnodes:
if isinstance(Bnode, functional.PointEvaluation):
# case: IntMom x PointEval
sd = ref_el.get_spatial_dimension()
pt_dict = {}
pt_old = Anode.get_point_dict()
for pt in pt_old:
pt_dict[pt+_first_point(Bnode)] = pt_old[pt]
# THE FOLLOWING IS PROBABLY CORRECT BUT UNTESTED
shp = (sd,)
nodes.append(functional.Functional(ref_el, shp, pt_dict, {}, "IntegralMoment"))
else:
> raise NotImplementedError("unsupported functional type")
E NotImplementedError: unsupported functional type
../firedrake_venv/src/fiat/FIAT/tensor_product.py:198: NotImplementedError
Check failure on line 57 in tests/regression/test_facet_split.py
github-actions / Firedrake real
test_facet_split.test_facet_split[jacobi-True]
NotImplementedError: unsupported functional type
Raw output
quadrilateral = True, pc_type = 'jacobi'
@pytest.mark.parametrize("quadrilateral", [True, False])
@pytest.mark.parametrize("pc_type", ["lu", "jacobi"])
def test_facet_split(quadrilateral, pc_type):
> assert run_facet_split(quadrilateral, pc_type) < 1E-10
tests/regression/test_facet_split.py:64:
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
tests/regression/test_facet_split.py:57: in run_facet_split
solve(a == L, uh, bcs=bcs, solver_parameters=parameters)
petsc4py/PETSc/Log.pyx:188: in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
???
petsc4py/PETSc/Log.pyx:189: in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
???
firedrake/adjoint_utils/solving.py:57: in wrapper
output = solve(*args, **kwargs)
firedrake/solving.py:140: in solve
_solve_varproblem(*args, **kwargs)
firedrake/solving.py:177: in _solve_varproblem
solver.solve()
petsc4py/PETSc/Log.pyx:188: in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
???
petsc4py/PETSc/Log.pyx:189: in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
???
firedrake/adjoint_utils/variational_solver.py:89: in wrapper
out = solve(self, **kwargs)
firedrake/variational_solver.py:307: in solve
dbc.apply(problem.u_restrict)
petsc4py/PETSc/Log.pyx:188: in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
???
petsc4py/PETSc/Log.pyx:189: in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
???
firedrake/adjoint_utils/dirichletbc.py:32: in wrapper
ret = apply(self, *args, **kwargs)
firedrake/bcs.py:451: in apply
r.assign(self.function_arg, subset=self.node_set)
../firedrake_venv/src/PyOP2/pyop2/utils.py:62: in __get__
obj.__dict__[self.__name__] = result = self.fget(obj)
firedrake/bcs.py:179: in node_set
return op2.Subset(self._function_space.node_set, self.nodes)
../firedrake_venv/src/PyOP2/pyop2/utils.py:62: in __get__
obj.__dict__[self.__name__] = result = self.fget(obj)
firedrake/bcs.py:169: in nodes
bcnodes1.append(hermite_stride(self._function_space.boundary_nodes(ss)))
firedrake/bcs.py:143: in hermite_stride
deriv_nodes = [i for i, node in enumerate(fe.fiat_equivalent.dual) if len(node.deriv_dict) != 0]
../firedrake_venv/src/tsfc/gem/utils.py:19: in __get__
obj.__dict__[self.__name__] = result = self.fget(obj)
../firedrake_venv/src/FInAT/finat/cube.py:65: in fiat_equivalent
return FIAT_FlattenedDimensions(self.product.fiat_equivalent)
../firedrake_venv/src/tsfc/gem/utils.py:19: in __get__
obj.__dict__[self.__name__] = result = self.fget(obj)
../firedrake_venv/src/FInAT/finat/tensor_product.py:82: in fiat_equivalent
return FIAT.TensorProductElement(A.fiat_equivalent, B.fiat_equivalent)
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
self = <FIAT.tensor_product.TensorProductElement object at 0x7fdfb38fee40>
A = <FIAT.fdm_element.FDMLagrange object at 0x7fdfb3469160>
B = <FIAT.fdm_element.FDMLagrange object at 0x7fdfb3469160>
def __init__(self, A, B):
# set up simple things
order = min(A.get_order(), B.get_order())
if A.get_formdegree() is None or B.get_formdegree() is None:
formdegree = None
else:
formdegree = A.get_formdegree() + B.get_formdegree()
# set up reference element
ref_el = TensorProductCell(A.get_reference_element(),
B.get_reference_element())
if A.mapping()[0] != "affine" and B.mapping()[0] == "affine":
mapping = A.mapping()[0]
elif B.mapping()[0] != "affine" and A.mapping()[0] == "affine":
mapping = B.mapping()[0]
elif A.mapping()[0] == "affine" and B.mapping()[0] == "affine":
mapping = "affine"
else:
raise ValueError("check tensor product mappings - at least one must be affine")
# set up entity_ids
Adofs = A.entity_dofs()
Bdofs = B.entity_dofs()
Bsdim = B.space_dimension()
entity_ids = {}
for curAdim in Adofs:
for curBdim in Bdofs:
entity_ids[(curAdim, curBdim)] = {}
dim_cur = 0
for entityA in Adofs[curAdim]:
for entityB in Bdofs[curBdim]:
entity_ids[(curAdim, curBdim)][dim_cur] = \
[x*Bsdim + y for x in Adofs[curAdim][entityA]
for y in Bdofs[curBdim][entityB]]
dim_cur += 1
# set up dual basis
Anodes = A.dual_basis()
Bnodes = B.dual_basis()
# build the dual set by inspecting the current dual
# sets item by item.
# Currently supported cases:
# PointEval x PointEval = PointEval [scalar x scalar = scalar]
# PointScaledNormalEval x PointEval = PointScaledNormalEval [vector x scalar = vector]
# ComponentPointEvaluation x PointEval [vector x scalar = vector]
nodes = []
for Anode in Anodes:
if isinstance(Anode, functional.PointEvaluation):
for Bnode in Bnodes:
if isinstance(Bnode, functional.PointEvaluation):
# case: PointEval x PointEval
# the PointEval functional just requires the
# coordinates. these are currently stored as
# the key of a one-item dictionary. we retrieve
# these by calling get_point_dict(), and
# use the concatenation to make a new PointEval
nodes.append(functional.PointEvaluation(ref_el, _first_point(Anode) + _first_point(Bnode)))
elif isinstance(Bnode, functional.IntegralMoment):
# dummy functional for product with integral moments
nodes.append(functional.Functional(None, None, None,
{}, "Undefined"))
elif isinstance(Bnode, functional.PointDerivative):
# dummy functional for product with point derivative
nodes.append(functional.Functional(None, None, None,
{}, "Undefined"))
else:
raise NotImplementedError("unsupported functional type")
elif isinstance(Anode, functional.PointScaledNormalEvaluation):
for Bnode in Bnodes:
if isinstance(Bnode, functional.PointEvaluation):
# case: PointScaledNormalEval x PointEval
# this could be wrong if the second shape
# has spatial dimension >1, since we are not
# explicitly scaling by facet size
if len(_first_point(Bnode)) > 1:
# TODO: support this case one day
raise NotImplementedError("PointScaledNormalEval x PointEval is not yet supported if the second shape has dimension > 1")
# We cannot make a new functional.PSNEval in
# the natural way, since it tries to compute
# the normal vector by itself.
# Instead, we create things manually, and
# call Functional() with these arguments
sd = ref_el.get_spatial_dimension()
# The pt_dict is a one-item dictionary containing
# the details of the functional.
# The key is the spatial coordinate, which
# is just a concatenation of the two parts.
# The value is a list of tuples, representing
# the normal vector (scaled by the volume of
# the facet) at that point.
# Each tuple looks like (foo, (i,)); the i'th
# component of the scaled normal is foo.
# The following line is only valid when the second
# shape has spatial dimension 1 (enforced above)
Apoint, Avalue = _first_point_pair(Anode)
pt_dict = {Apoint + _first_point(Bnode): Avalue + [(0.0, (len(Apoint),))]}
# The following line should be used in the
# general case
# pt_dict = {Anode.get_point_dict().keys()[0] + Bnode.get_point_dict().keys()[0]: Anode.get_point_dict().values()[0] + [(0.0, (ii,)) for ii in range(len(Anode.get_point_dict().keys()[0]), len(Anode.get_point_dict().keys()[0]) + len(Bnode.get_point_dict().keys()[0]))]}
# THE FOLLOWING IS PROBABLY CORRECT BUT UNTESTED
shp = (sd,)
nodes.append(functional.Functional(ref_el, shp, pt_dict, {}, "PointScaledNormalEval"))
else:
raise NotImplementedError("unsupported functional type")
elif isinstance(Anode, functional.PointEdgeTangentEvaluation):
for Bnode in Bnodes:
if isinstance(Bnode, functional.PointEvaluation):
# case: PointEdgeTangentEval x PointEval
# this is very similar to the case above, so comments omitted
if len(_first_point(Bnode)) > 1:
raise NotImplementedError("PointEdgeTangentEval x PointEval is not yet supported if the second shape has dimension > 1")
sd = ref_el.get_spatial_dimension()
Apoint, Avalue = _first_point_pair(Anode)
pt_dict = {Apoint + _first_point(Bnode): Avalue + [(0.0, (len(Apoint),))]}
# THE FOLLOWING IS PROBABLY CORRECT BUT UNTESTED
shp = (sd,)
nodes.append(functional.Functional(ref_el, shp, pt_dict, {}, "PointEdgeTangent"))
else:
raise NotImplementedError("unsupported functional type")
elif isinstance(Anode, functional.ComponentPointEvaluation):
for Bnode in Bnodes:
if isinstance(Bnode, functional.PointEvaluation):
# case: ComponentPointEval x PointEval
# the CptPointEval functional requires the component
# and the coordinates. very similar to PE x PE case.
sd = ref_el.get_spatial_dimension()
nodes.append(functional.ComponentPointEvaluation(ref_el, Anode.comp, (sd,), _first_point(Anode) + _first_point(Bnode)))
else:
raise NotImplementedError("unsupported functional type")
elif isinstance(Anode, functional.FrobeniusIntegralMoment):
for Bnode in Bnodes:
if isinstance(Bnode, functional.PointEvaluation):
# case: FroIntMom x PointEval
sd = ref_el.get_spatial_dimension()
pt_dict = {}
pt_old = Anode.get_point_dict()
for pt in pt_old:
pt_dict[pt+_first_point(Bnode)] = pt_old[pt] + [(0.0, sd-1)]
# THE FOLLOWING IS PROBABLY CORRECT BUT UNTESTED
shp = (sd,)
nodes.append(functional.Functional(ref_el, shp, pt_dict, {}, "FrobeniusIntegralMoment"))
else:
raise NotImplementedError("unsupported functional type")
elif isinstance(Anode, functional.IntegralMoment):
for Bnode in Bnodes:
if isinstance(Bnode, functional.PointEvaluation):
# case: IntMom x PointEval
sd = ref_el.get_spatial_dimension()
pt_dict = {}
pt_old = Anode.get_point_dict()
for pt in pt_old:
pt_dict[pt+_first_point(Bnode)] = pt_old[pt]
# THE FOLLOWING IS PROBABLY CORRECT BUT UNTESTED
shp = (sd,)
nodes.append(functional.Functional(ref_el, shp, pt_dict, {}, "IntegralMoment"))
else:
> raise NotImplementedError("unsupported functional type")
E NotImplementedError: unsupported functional type
../firedrake_venv/src/fiat/FIAT/tensor_product.py:198: NotImplementedError
Check failure on line 1 in tests/regression/test_facet_split.py
github-actions / Firedrake real
test_facet_split.test_facet_split_parallel[lu]
subprocess.CalledProcessError: Command '['mpiexec', '-n', '1', '-genv', '_PYTEST_MPI_CHILD_PROCESS', '1', '/__w/firedrake/firedrake_venv/bin/python', '-m', 'pytest', '--runxfail', '-s', '-q', '/__w/firedrake/firedrake/tests/regression/test_facet_split.py::test_facet_split_parallel[lu]', ':', '-n', '2', '/__w/firedrake/firedrake_venv/bin/python', '-m', 'pytest', '--runxfail', '-s', '-q', '/__w/firedrake/firedrake/tests/regression/test_facet_split.py::test_facet_split_parallel[lu]', '--tb=no', '--no-summary', '--no-header', '--disable-warnings', '--show-capture=no']' returned non-zero exit status 1.
Raw output
args = (), kwargs = {'pc_type': 'lu'}
def parallel_callback(*args, **kwargs):
> subprocess.run(cmd, check=True)
../firedrake_venv/src/pytest-mpi/pytest_mpi.py:210:
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
input = None, capture_output = False, timeout = None, check = True
popenargs = (['mpiexec', '-n', '1', '-genv', '_PYTEST_MPI_CHILD_PROCESS', '1', ...],)
kwargs = {}
process = <Popen: returncode: 1 args: ['mpiexec', '-n', '1', '-genv', '_PYTEST_MPI_CHI...>
stdout = None, stderr = None, retcode = 1
def run(*popenargs,
input=None, capture_output=False, timeout=None, check=False, **kwargs):
"""Run command with arguments and return a CompletedProcess instance.
The returned instance will have attributes args, returncode, stdout and
stderr. By default, stdout and stderr are not captured, and those attributes
will be None. Pass stdout=PIPE and/or stderr=PIPE in order to capture them,
or pass capture_output=True to capture both.
If check is True and the exit code was non-zero, it raises a
CalledProcessError. The CalledProcessError object will have the return code
in the returncode attribute, and output & stderr attributes if those streams
were captured.
If timeout is given, and the process takes too long, a TimeoutExpired
exception will be raised.
There is an optional argument "input", allowing you to
pass bytes or a string to the subprocess's stdin. If you use this argument
you may not also use the Popen constructor's "stdin" argument, as
it will be used internally.
By default, all communication is in bytes, and therefore any "input" should
be bytes, and the stdout and stderr will be bytes. If in text mode, any
"input" should be a string, and stdout and stderr will be strings decoded
according to locale encoding, or by "encoding" if set. Text mode is
triggered by setting any of text, encoding, errors or universal_newlines.
The other arguments are the same as for the Popen constructor.
"""
if input is not None:
if kwargs.get('stdin') is not None:
raise ValueError('stdin and input arguments may not both be used.')
kwargs['stdin'] = PIPE
if capture_output:
if kwargs.get('stdout') is not None or kwargs.get('stderr') is not None:
raise ValueError('stdout and stderr arguments may not be used '
'with capture_output.')
kwargs['stdout'] = PIPE
kwargs['stderr'] = PIPE
with Popen(*popenargs, **kwargs) as process:
try:
stdout, stderr = process.communicate(input, timeout=timeout)
except TimeoutExpired as exc:
process.kill()
if _mswindows:
# Windows accumulates the output in a single blocking
# read() call run on child threads, with the timeout
# being done in a join() on those threads. communicate()
# _after_ kill() is required to collect that and add it
# to the exception.
exc.stdout, exc.stderr = process.communicate()
else:
# POSIX _communicate already populated the output so
# far into the TimeoutExpired exception.
process.wait()
raise
except: # Including KeyboardInterrupt, communicate handled that.
process.kill()
# We don't call process.wait() as .__exit__ does that for us.
raise
retcode = process.poll()
if check and retcode:
> raise CalledProcessError(retcode, process.args,
output=stdout, stderr=stderr)
E subprocess.CalledProcessError: Command '['mpiexec', '-n', '1', '-genv', '_PYTEST_MPI_CHILD_PROCESS', '1', '/__w/firedrake/firedrake_venv/bin/python', '-m', 'pytest', '--runxfail', '-s', '-q', '/__w/firedrake/firedrake/tests/regression/test_facet_split.py::test_facet_split_parallel[lu]', ':', '-n', '2', '/__w/firedrake/firedrake_venv/bin/python', '-m', 'pytest', '--runxfail', '-s', '-q', '/__w/firedrake/firedrake/tests/regression/test_facet_split.py::test_facet_split_parallel[lu]', '--tb=no', '--no-summary', '--no-header', '--disable-warnings', '--show-capture=no']' returned non-zero exit status 1.
/usr/lib/python3.12/subprocess.py:571: CalledProcessError
Check failure on line 1 in tests/regression/test_facet_split.py
github-actions / Firedrake real
test_facet_split.test_facet_split_parallel[jacobi]
subprocess.CalledProcessError: Command '['mpiexec', '-n', '1', '-genv', '_PYTEST_MPI_CHILD_PROCESS', '1', '/__w/firedrake/firedrake_venv/bin/python', '-m', 'pytest', '--runxfail', '-s', '-q', '/__w/firedrake/firedrake/tests/regression/test_facet_split.py::test_facet_split_parallel[jacobi]', ':', '-n', '2', '/__w/firedrake/firedrake_venv/bin/python', '-m', 'pytest', '--runxfail', '-s', '-q', '/__w/firedrake/firedrake/tests/regression/test_facet_split.py::test_facet_split_parallel[jacobi]', '--tb=no', '--no-summary', '--no-header', '--disable-warnings', '--show-capture=no']' returned non-zero exit status 1.
Raw output
args = (), kwargs = {'pc_type': 'jacobi'}
def parallel_callback(*args, **kwargs):
> subprocess.run(cmd, check=True)
../firedrake_venv/src/pytest-mpi/pytest_mpi.py:210:
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
input = None, capture_output = False, timeout = None, check = True
popenargs = (['mpiexec', '-n', '1', '-genv', '_PYTEST_MPI_CHILD_PROCESS', '1', ...],)
kwargs = {}
process = <Popen: returncode: 1 args: ['mpiexec', '-n', '1', '-genv', '_PYTEST_MPI_CHI...>
stdout = None, stderr = None, retcode = 1
def run(*popenargs,
input=None, capture_output=False, timeout=None, check=False, **kwargs):
"""Run command with arguments and return a CompletedProcess instance.
The returned instance will have attributes args, returncode, stdout and
stderr. By default, stdout and stderr are not captured, and those attributes
will be None. Pass stdout=PIPE and/or stderr=PIPE in order to capture them,
or pass capture_output=True to capture both.
If check is True and the exit code was non-zero, it raises a
CalledProcessError. The CalledProcessError object will have the return code
in the returncode attribute, and output & stderr attributes if those streams
were captured.
If timeout is given, and the process takes too long, a TimeoutExpired
exception will be raised.
There is an optional argument "input", allowing you to
pass bytes or a string to the subprocess's stdin. If you use this argument
you may not also use the Popen constructor's "stdin" argument, as
it will be used internally.
By default, all communication is in bytes, and therefore any "input" should
be bytes, and the stdout and stderr will be bytes. If in text mode, any
"input" should be a string, and stdout and stderr will be strings decoded
according to locale encoding, or by "encoding" if set. Text mode is
triggered by setting any of text, encoding, errors or universal_newlines.
The other arguments are the same as for the Popen constructor.
"""
if input is not None:
if kwargs.get('stdin') is not None:
raise ValueError('stdin and input arguments may not both be used.')
kwargs['stdin'] = PIPE
if capture_output:
if kwargs.get('stdout') is not None or kwargs.get('stderr') is not None:
raise ValueError('stdout and stderr arguments may not be used '
'with capture_output.')
kwargs['stdout'] = PIPE
kwargs['stderr'] = PIPE
with Popen(*popenargs, **kwargs) as process:
try:
stdout, stderr = process.communicate(input, timeout=timeout)
except TimeoutExpired as exc:
process.kill()
if _mswindows:
# Windows accumulates the output in a single blocking
# read() call run on child threads, with the timeout
# being done in a join() on those threads. communicate()
# _after_ kill() is required to collect that and add it
# to the exception.
exc.stdout, exc.stderr = process.communicate()
else:
# POSIX _communicate already populated the output so
# far into the TimeoutExpired exception.
process.wait()
raise
except: # Including KeyboardInterrupt, communicate handled that.
process.kill()
# We don't call process.wait() as .__exit__ does that for us.
raise
retcode = process.poll()
if check and retcode:
> raise CalledProcessError(retcode, process.args,
output=stdout, stderr=stderr)
E subprocess.CalledProcessError: Command '['mpiexec', '-n', '1', '-genv', '_PYTEST_MPI_CHILD_PROCESS', '1', '/__w/firedrake/firedrake_venv/bin/python', '-m', 'pytest', '--runxfail', '-s', '-q', '/__w/firedrake/firedrake/tests/regression/test_facet_split.py::test_facet_split_parallel[jacobi]', ':', '-n', '2', '/__w/firedrake/firedrake_venv/bin/python', '-m', 'pytest', '--runxfail', '-s', '-q', '/__w/firedrake/firedrake/tests/regression/test_facet_split.py::test_facet_split_parallel[jacobi]', '--tb=no', '--no-summary', '--no-header', '--disable-warnings', '--show-capture=no']' returned non-zero exit status 1.
/usr/lib/python3.12/subprocess.py:571: CalledProcessError
Check failure on line 114 in tests/regression/test_fdm.py
github-actions / Firedrake real
test_fdm.test_p_independence_hgrad[Rectangle-spectral]
NotImplementedError: unsupported functional type
Raw output
mesh = Mesh(VectorElement(FiniteElement('Q', quadrilateral, 1), dim=2), 3186)
variant = None
@pytest.mark.skipcomplex
def test_p_independence_hgrad(mesh, variant):
family = "Lagrange"
expected = [16, 12] if mesh.topological_dimension() == 3 else [9, 7]
solvers = [fdmstar] if variant is None else [fdmstar, facetstar]
for degree in range(3, 6):
V = FunctionSpace(mesh, family, degree, variant=variant)
problem = build_riesz_map(V, grad)
for sp, expected_it in zip(solvers, expected):
> assert solve_riesz_map(problem, sp) <= expected_it
tests/regression/test_fdm.py:147:
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
tests/regression/test_fdm.py:114: in solve_riesz_map
solver.solve()
petsc4py/PETSc/Log.pyx:188: in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
???
petsc4py/PETSc/Log.pyx:189: in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
???
firedrake/adjoint_utils/variational_solver.py:89: in wrapper
out = solve(self, **kwargs)
firedrake/variational_solver.py:324: in solve
self.snes.solve(None, work)
petsc4py/PETSc/SNES.pyx:1601: in petsc4py.PETSc.SNES.solve
???
petsc4py/PETSc/PETSc.pyx:323: in petsc4py.PETSc.PetscPythonErrorHandler
???
petsc4py/PETSc/PETSc.pyx:323: in petsc4py.PETSc.PetscPythonErrorHandler
???
petsc4py/PETSc/PETSc.pyx:323: in petsc4py.PETSc.PetscPythonErrorHandler
???
petsc4py/PETSc/PETSc.pyx:323: in petsc4py.PETSc.PetscPythonErrorHandler
???
petsc4py/PETSc/PETSc.pyx:323: in petsc4py.PETSc.PetscPythonErrorHandler
???
petsc4py/PETSc/PETSc.pyx:323: in petsc4py.PETSc.PetscPythonErrorHandler
???
petsc4py/PETSc/libpetsc4py.pyx:1324: in petsc4py.PETSc.PCSetUp_Python
???
firedrake/preconditioners/base.py:129: in setUp
super().setUp(pc)
firedrake/preconditioners/base.py:45: in setUp
self.initialize(pc)
firedrake/preconditioners/pmg.py:157: in initialize
ppc.setUp()
petsc4py/PETSc/PC.pyx:563: in petsc4py.PETSc.PC.setUp
???
petsc4py/PETSc/PETSc.pyx:323: in petsc4py.PETSc.PetscPythonErrorHandler
???
petsc4py/PETSc/PETSc.pyx:323: in petsc4py.PETSc.PetscPythonErrorHandler
???
petsc4py/PETSc/PETSc.pyx:323: in petsc4py.PETSc.PetscPythonErrorHandler
???
petsc4py/PETSc/PETSc.pyx:323: in petsc4py.PETSc.PetscPythonErrorHandler
???
petsc4py/PETSc/PETSc.pyx:323: in petsc4py.PETSc.PetscPythonErrorHandler
???
petsc4py/PETSc/PETSc.pyx:323: in petsc4py.PETSc.PetscPythonErrorHandler
???
petsc4py/PETSc/PETSc.pyx:323: in petsc4py.PETSc.PetscPythonErrorHandler
???
petsc4py/PETSc/PETSc.pyx:323: in petsc4py.PETSc.PetscPythonErrorHandler
???
petsc4py/PETSc/libpetsc4py.pyx:1324: in petsc4py.PETSc.PCSetUp_Python
???
firedrake/preconditioners/base.py:129: in setUp
super().setUp(pc)
firedrake/preconditioners/base.py:45: in setUp
self.initialize(pc)
petsc4py/PETSc/Log.pyx:188: in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
???
petsc4py/PETSc/Log.pyx:189: in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
???
firedrake/preconditioners/fdm.py:178: in initialize
Amat, Pmat, self.assembly_callables = self.allocate_matrix(Amat, V_fdm, J_fdm, bcs_fdm, fcp,
petsc4py/PETSc/Log.pyx:188: in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
???
petsc4py/PETSc/Log.pyx:189: in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
???
firedrake/preconditioners/fdm.py:235: in allocate_matrix
self.lgmaps = {Vsub: Vsub.local_to_global_map([bc for bc in bcs if bc.function_space() == Vsub]) for Vsub in V}
petsc4py/PETSc/Log.pyx:188: in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
???
petsc4py/PETSc/Log.pyx:189: in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
???
firedrake/functionspaceimpl.py:840: in local_to_global_map
nodes.append(bc.nodes)
../firedrake_venv/src/PyOP2/pyop2/utils.py:62: in __get__
obj.__dict__[self.__name__] = result = self.fget(obj)
firedrake/bcs.py:169: in nodes
bcnodes1.append(hermite_stride(self._function_space.boundary_nodes(ss)))
firedrake/bcs.py:143: in hermite_stride
deriv_nodes = [i for i, node in enumerate(fe.fiat_equivalent.dual) if len(node.deriv_dict) != 0]
../firedrake_venv/src/tsfc/gem/utils.py:19: in __get__
obj.__dict__[self.__name__] = result = self.fget(obj)
../firedrake_venv/src/FInAT/finat/cube.py:65: in fiat_equivalent
return FIAT_FlattenedDimensions(self.product.fiat_equivalent)
../firedrake_venv/src/tsfc/gem/utils.py:19: in __get__
obj.__dict__[self.__name__] = result = self.fget(obj)
../firedrake_venv/src/FInAT/finat/tensor_product.py:82: in fiat_equivalent
return FIAT.TensorProductElement(A.fiat_equivalent, B.fiat_equivalent)
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
self = <FIAT.tensor_product.TensorProductElement object at 0x7fdfb2ea2e10>
A = <FIAT.fdm_element.FDMLagrange object at 0x7fdfb3469160>
B = <FIAT.fdm_element.FDMLagrange object at 0x7fdfb3469160>
def __init__(self, A, B):
# set up simple things
order = min(A.get_order(), B.get_order())
if A.get_formdegree() is None or B.get_formdegree() is None:
formdegree = None
else:
formdegree = A.get_formdegree() + B.get_formdegree()
# set up reference element
ref_el = TensorProductCell(A.get_reference_element(),
B.get_reference_element())
if A.mapping()[0] != "affine" and B.mapping()[0] == "affine":
mapping = A.mapping()[0]
elif B.mapping()[0] != "affine" and A.mapping()[0] == "affine":
mapping = B.mapping()[0]
elif A.mapping()[0] == "affine" and B.mapping()[0] == "affine":
mapping = "affine"
else:
raise ValueError("check tensor product mappings - at least one must be affine")
# set up entity_ids
Adofs = A.entity_dofs()
Bdofs = B.entity_dofs()
Bsdim = B.space_dimension()
entity_ids = {}
for curAdim in Adofs:
for curBdim in Bdofs:
entity_ids[(curAdim, curBdim)] = {}
dim_cur = 0
for entityA in Adofs[curAdim]:
for entityB in Bdofs[curBdim]:
entity_ids[(curAdim, curBdim)][dim_cur] = \
[x*Bsdim + y for x in Adofs[curAdim][entityA]
for y in Bdofs[curBdim][entityB]]
dim_cur += 1
# set up dual basis
Anodes = A.dual_basis()
Bnodes = B.dual_basis()
# build the dual set by inspecting the current dual
# sets item by item.
# Currently supported cases:
# PointEval x PointEval = PointEval [scalar x scalar = scalar]
# PointScaledNormalEval x PointEval = PointScaledNormalEval [vector x scalar = vector]
# ComponentPointEvaluation x PointEval [vector x scalar = vector]
nodes = []
for Anode in Anodes:
if isinstance(Anode, functional.PointEvaluation):
for Bnode in Bnodes:
if isinstance(Bnode, functional.PointEvaluation):
# case: PointEval x PointEval
# the PointEval functional just requires the
# coordinates. these are currently stored as
# the key of a one-item dictionary. we retrieve
# these by calling get_point_dict(), and
# use the concatenation to make a new PointEval
nodes.append(functional.PointEvaluation(ref_el, _first_point(Anode) + _first_point(Bnode)))
elif isinstance(Bnode, functional.IntegralMoment):
# dummy functional for product with integral moments
nodes.append(functional.Functional(None, None, None,
{}, "Undefined"))
elif isinstance(Bnode, functional.PointDerivative):
# dummy functional for product with point derivative
nodes.append(functional.Functional(None, None, None,
{}, "Undefined"))
else:
raise NotImplementedError("unsupported functional type")
elif isinstance(Anode, functional.PointScaledNormalEvaluation):
for Bnode in Bnodes:
if isinstance(Bnode, functional.PointEvaluation):
# case: PointScaledNormalEval x PointEval
# this could be wrong if the second shape
# has spatial dimension >1, since we are not
# explicitly scaling by facet size
if len(_first_point(Bnode)) > 1:
# TODO: support this case one day
raise NotImplementedError("PointScaledNormalEval x PointEval is not yet supported if the second shape has dimension > 1")
# We cannot make a new functional.PSNEval in
# the natural way, since it tries to compute
# the normal vector by itself.
# Instead, we create things manually, and
# call Functional() with these arguments
sd = ref_el.get_spatial_dimension()
# The pt_dict is a one-item dictionary containing
# the details of the functional.
# The key is the spatial coordinate, which
# is just a concatenation of the two parts.
# The value is a list of tuples, representing
# the normal vector (scaled by the volume of
# the facet) at that point.
# Each tuple looks like (foo, (i,)); the i'th
# component of the scaled normal is foo.
# The following line is only valid when the second
# shape has spatial dimension 1 (enforced above)
Apoint, Avalue = _first_point_pair(Anode)
pt_dict = {Apoint + _first_point(Bnode): Avalue + [(0.0, (len(Apoint),))]}
# The following line should be used in the
# general case
# pt_dict = {Anode.get_point_dict().keys()[0] + Bnode.get_point_dict().keys()[0]: Anode.get_point_dict().values()[0] + [(0.0, (ii,)) for ii in range(len(Anode.get_point_dict().keys()[0]), len(Anode.get_point_dict().keys()[0]) + len(Bnode.get_point_dict().keys()[0]))]}
# THE FOLLOWING IS PROBABLY CORRECT BUT UNTESTED
shp = (sd,)
nodes.append(functional.Functional(ref_el, shp, pt_dict, {}, "PointScaledNormalEval"))
else:
raise NotImplementedError("unsupported functional type")
elif isinstance(Anode, functional.PointEdgeTangentEvaluation):
for Bnode in Bnodes:
if isinstance(Bnode, functional.PointEvaluation):
# case: PointEdgeTangentEval x PointEval
# this is very similar to the case above, so comments omitted
if len(_first_point(Bnode)) > 1:
raise NotImplementedError("PointEdgeTangentEval x PointEval is not yet supported if the second shape has dimension > 1")
sd = ref_el.get_spatial_dimension()
Apoint, Avalue = _first_point_pair(Anode)
pt_dict = {Apoint + _first_point(Bnode): Avalue + [(0.0, (len(Apoint),))]}
# THE FOLLOWING IS PROBABLY CORRECT BUT UNTESTED
shp = (sd,)
nodes.append(functional.Functional(ref_el, shp, pt_dict, {}, "PointEdgeTangent"))
else:
raise NotImplementedError("unsupported functional type")
elif isinstance(Anode, functional.ComponentPointEvaluation):
for Bnode in Bnodes:
if isinstance(Bnode, functional.PointEvaluation):
# case: ComponentPointEval x PointEval
# the CptPointEval functional requires the component
# and the coordinates. very similar to PE x PE case.
sd = ref_el.get_spatial_dimension()
nodes.append(functional.ComponentPointEvaluation(ref_el, Anode.comp, (sd,), _first_point(Anode) + _first_point(Bnode)))
else:
raise NotImplementedError("unsupported functional type")
elif isinstance(Anode, functional.FrobeniusIntegralMoment):
for Bnode in Bnodes:
if isinstance(Bnode, functional.PointEvaluation):
# case: FroIntMom x PointEval
sd = ref_el.get_spatial_dimension()
pt_dict = {}
pt_old = Anode.get_point_dict()
for pt in pt_old:
pt_dict[pt+_first_point(Bnode)] = pt_old[pt] + [(0.0, sd-1)]
# THE FOLLOWING IS PROBABLY CORRECT BUT UNTESTED
shp = (sd,)
nodes.append(functional.Functional(ref_el, shp, pt_dict, {}, "FrobeniusIntegralMoment"))
else:
raise NotImplementedError("unsupported functional type")
elif isinstance(Anode, functional.IntegralMoment):
for Bnode in Bnodes:
if isinstance(Bnode, functional.PointEvaluation):
# case: IntMom x PointEval
sd = ref_el.get_spatial_dimension()
pt_dict = {}
pt_old = Anode.get_point_dict()
for pt in pt_old:
pt_dict[pt+_first_point(Bnode)] = pt_old[pt]
# THE FOLLOWING IS PROBABLY CORRECT BUT UNTESTED
shp = (sd,)
nodes.append(functional.Functional(ref_el, shp, pt_dict, {}, "IntegralMoment"))
else:
> raise NotImplementedError("unsupported functional type")
E NotImplementedError: unsupported functional type
../firedrake_venv/src/fiat/FIAT/tensor_product.py:198: NotImplementedError
Check failure on line 114 in tests/regression/test_fdm.py
github-actions / Firedrake real
test_fdm.test_p_independence_hgrad[Rectangle-fdm]
NotImplementedError: unsupported functional type
Raw output
mesh = Mesh(VectorElement(FiniteElement('Q', quadrilateral, 1), dim=2), 3212)
variant = 'fdm'
@pytest.mark.skipcomplex
def test_p_independence_hgrad(mesh, variant):
family = "Lagrange"
expected = [16, 12] if mesh.topological_dimension() == 3 else [9, 7]
solvers = [fdmstar] if variant is None else [fdmstar, facetstar]
for degree in range(3, 6):
V = FunctionSpace(mesh, family, degree, variant=variant)
problem = build_riesz_map(V, grad)
for sp, expected_it in zip(solvers, expected):
> assert solve_riesz_map(problem, sp) <= expected_it
tests/regression/test_fdm.py:147:
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
tests/regression/test_fdm.py:114: in solve_riesz_map
solver.solve()
petsc4py/PETSc/Log.pyx:188: in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
???
petsc4py/PETSc/Log.pyx:189: in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
???
firedrake/adjoint_utils/variational_solver.py:89: in wrapper
out = solve(self, **kwargs)
firedrake/variational_solver.py:307: in solve
dbc.apply(problem.u_restrict)
petsc4py/PETSc/Log.pyx:188: in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
???
petsc4py/PETSc/Log.pyx:189: in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
???
firedrake/adjoint_utils/dirichletbc.py:32: in wrapper
ret = apply(self, *args, **kwargs)
firedrake/bcs.py:451: in apply
r.assign(self.function_arg, subset=self.node_set)
../firedrake_venv/src/PyOP2/pyop2/utils.py:62: in __get__
obj.__dict__[self.__name__] = result = self.fget(obj)
firedrake/bcs.py:179: in node_set
return op2.Subset(self._function_space.node_set, self.nodes)
../firedrake_venv/src/PyOP2/pyop2/utils.py:62: in __get__
obj.__dict__[self.__name__] = result = self.fget(obj)
firedrake/bcs.py:169: in nodes
bcnodes1.append(hermite_stride(self._function_space.boundary_nodes(ss)))
firedrake/bcs.py:143: in hermite_stride
deriv_nodes = [i for i, node in enumerate(fe.fiat_equivalent.dual) if len(node.deriv_dict) != 0]
../firedrake_venv/src/tsfc/gem/utils.py:19: in __get__
obj.__dict__[self.__name__] = result = self.fget(obj)
../firedrake_venv/src/FInAT/finat/cube.py:65: in fiat_equivalent
return FIAT_FlattenedDimensions(self.product.fiat_equivalent)
../firedrake_venv/src/tsfc/gem/utils.py:19: in __get__
obj.__dict__[self.__name__] = result = self.fget(obj)
../firedrake_venv/src/FInAT/finat/tensor_product.py:82: in fiat_equivalent
return FIAT.TensorProductElement(A.fiat_equivalent, B.fiat_equivalent)
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
self = <FIAT.tensor_product.TensorProductElement object at 0x7fdfb3d3acf0>
A = <FIAT.fdm_element.FDMLagrange object at 0x7fdfb3469160>
B = <FIAT.fdm_element.FDMLagrange object at 0x7fdfb3469160>
def __init__(self, A, B):
# set up simple things
order = min(A.get_order(), B.get_order())
if A.get_formdegree() is None or B.get_formdegree() is None:
formdegree = None
else:
formdegree = A.get_formdegree() + B.get_formdegree()
# set up reference element
ref_el = TensorProductCell(A.get_reference_element(),
B.get_reference_element())
if A.mapping()[0] != "affine" and B.mapping()[0] == "affine":
mapping = A.mapping()[0]
elif B.mapping()[0] != "affine" and A.mapping()[0] == "affine":
mapping = B.mapping()[0]
elif A.mapping()[0] == "affine" and B.mapping()[0] == "affine":
mapping = "affine"
else:
raise ValueError("check tensor product mappings - at least one must be affine")
# set up entity_ids
Adofs = A.entity_dofs()
Bdofs = B.entity_dofs()
Bsdim = B.space_dimension()
entity_ids = {}
for curAdim in Adofs:
for curBdim in Bdofs:
entity_ids[(curAdim, curBdim)] = {}
dim_cur = 0
for entityA in Adofs[curAdim]:
for entityB in Bdofs[curBdim]:
entity_ids[(curAdim, curBdim)][dim_cur] = \
[x*Bsdim + y for x in Adofs[curAdim][entityA]
for y in Bdofs[curBdim][entityB]]
dim_cur += 1
# set up dual basis
Anodes = A.dual_basis()
Bnodes = B.dual_basis()
# build the dual set by inspecting the current dual
# sets item by item.
# Currently supported cases:
# PointEval x PointEval = PointEval [scalar x scalar = scalar]
# PointScaledNormalEval x PointEval = PointScaledNormalEval [vector x scalar = vector]
# ComponentPointEvaluation x PointEval [vector x scalar = vector]
nodes = []
for Anode in Anodes:
if isinstance(Anode, functional.PointEvaluation):
for Bnode in Bnodes:
if isinstance(Bnode, functional.PointEvaluation):
# case: PointEval x PointEval
# the PointEval functional just requires the
# coordinates. these are currently stored as
# the key of a one-item dictionary. we retrieve
# these by calling get_point_dict(), and
# use the concatenation to make a new PointEval
nodes.append(functional.PointEvaluation(ref_el, _first_point(Anode) + _first_point(Bnode)))
elif isinstance(Bnode, functional.IntegralMoment):
# dummy functional for product with integral moments
nodes.append(functional.Functional(None, None, None,
{}, "Undefined"))
elif isinstance(Bnode, functional.PointDerivative):
# dummy functional for product with point derivative
nodes.append(functional.Functional(None, None, None,
{}, "Undefined"))
else:
raise NotImplementedError("unsupported functional type")
elif isinstance(Anode, functional.PointScaledNormalEvaluation):
for Bnode in Bnodes:
if isinstance(Bnode, functional.PointEvaluation):
# case: PointScaledNormalEval x PointEval
# this could be wrong if the second shape
# has spatial dimension >1, since we are not
# explicitly scaling by facet size
if len(_first_point(Bnode)) > 1:
# TODO: support this case one day
raise NotImplementedError("PointScaledNormalEval x PointEval is not yet supported if the second shape has dimension > 1")
# We cannot make a new functional.PSNEval in
# the natural way, since it tries to compute
# the normal vector by itself.
# Instead, we create things manually, and
# call Functional() with these arguments
sd = ref_el.get_spatial_dimension()
# The pt_dict is a one-item dictionary containing
# the details of the functional.
# The key is the spatial coordinate, which
# is just a concatenation of the two parts.
# The value is a list of tuples, representing
# the normal vector (scaled by the volume of
# the facet) at that point.
# Each tuple looks like (foo, (i,)); the i'th
# component of the scaled normal is foo.
# The following line is only valid when the second
# shape has spatial dimension 1 (enforced above)
Apoint, Avalue = _first_point_pair(Anode)
pt_dict = {Apoint + _first_point(Bnode): Avalue + [(0.0, (len(Apoint),))]}
# The following line should be used in the
# general case
# pt_dict = {Anode.get_point_dict().keys()[0] + Bnode.get_point_dict().keys()[0]: Anode.get_point_dict().values()[0] + [(0.0, (ii,)) for ii in range(len(Anode.get_point_dict().keys()[0]), len(Anode.get_point_dict().keys()[0]) + len(Bnode.get_point_dict().keys()[0]))]}
# THE FOLLOWING IS PROBABLY CORRECT BUT UNTESTED
shp = (sd,)
nodes.append(functional.Functional(ref_el, shp, pt_dict, {}, "PointScaledNormalEval"))
else:
raise NotImplementedError("unsupported functional type")
elif isinstance(Anode, functional.PointEdgeTangentEvaluation):
for Bnode in Bnodes:
if isinstance(Bnode, functional.PointEvaluation):
# case: PointEdgeTangentEval x PointEval
# this is very similar to the case above, so comments omitted
if len(_first_point(Bnode)) > 1:
raise NotImplementedError("PointEdgeTangentEval x PointEval is not yet supported if the second shape has dimension > 1")
sd = ref_el.get_spatial_dimension()
Apoint, Avalue = _first_point_pair(Anode)
pt_dict = {Apoint + _first_point(Bnode): Avalue + [(0.0, (len(Apoint),))]}
# THE FOLLOWING IS PROBABLY CORRECT BUT UNTESTED
shp = (sd,)
nodes.append(functional.Functional(ref_el, shp, pt_dict, {}, "PointEdgeTangent"))
else:
raise NotImplementedError("unsupported functional type")
elif isinstance(Anode, functional.ComponentPointEvaluation):
for Bnode in Bnodes:
if isinstance(Bnode, functional.PointEvaluation):
# case: ComponentPointEval x PointEval
# the CptPointEval functional requires the component
# and the coordinates. very similar to PE x PE case.
sd = ref_el.get_spatial_dimension()
nodes.append(functional.ComponentPointEvaluation(ref_el, Anode.comp, (sd,), _first_point(Anode) + _first_point(Bnode)))
else:
raise NotImplementedError("unsupported functional type")
elif isinstance(Anode, functional.FrobeniusIntegralMoment):
for Bnode in Bnodes:
if isinstance(Bnode, functional.PointEvaluation):
# case: FroIntMom x PointEval
sd = ref_el.get_spatial_dimension()
pt_dict = {}
pt_old = Anode.get_point_dict()
for pt in pt_old:
pt_dict[pt+_first_point(Bnode)] = pt_old[pt] + [(0.0, sd-1)]
# THE FOLLOWING IS PROBABLY CORRECT BUT UNTESTED
shp = (sd,)
nodes.append(functional.Functional(ref_el, shp, pt_dict, {}, "FrobeniusIntegralMoment"))
else:
raise NotImplementedError("unsupported functional type")
elif isinstance(Anode, functional.IntegralMoment):
for Bnode in Bnodes:
if isinstance(Bnode, functional.PointEvaluation):
# case: IntMom x PointEval
sd = ref_el.get_spatial_dimension()
pt_dict = {}
pt_old = Anode.get_point_dict()
for pt in pt_old:
pt_dict[pt+_first_point(Bnode)] = pt_old[pt]
# THE FOLLOWING IS PROBABLY CORRECT BUT UNTESTED
shp = (sd,)
nodes.append(functional.Functional(ref_el, shp, pt_dict, {}, "IntegralMoment"))
else:
> raise NotImplementedError("unsupported functional type")
E NotImplementedError: unsupported functional type
../firedrake_venv/src/fiat/FIAT/tensor_product.py:198: NotImplementedError
Check failure on line 114 in tests/regression/test_fdm.py
github-actions / Firedrake real
test_fdm.test_p_independence_hgrad[Box-spectral]
NotImplementedError: unsupported functional type
Raw output
mesh = Mesh(VectorElement(TensorProductElement(FiniteElement('Q', quadrilateral, 1), FiniteElement('Lagrange', interval, 1), cell=TensorProductCell(quadrilateral, interval)), dim=3), 3228)
variant = None
@pytest.mark.skipcomplex
def test_p_independence_hgrad(mesh, variant):
family = "Lagrange"
expected = [16, 12] if mesh.topological_dimension() == 3 else [9, 7]
solvers = [fdmstar] if variant is None else [fdmstar, facetstar]
for degree in range(3, 6):
V = FunctionSpace(mesh, family, degree, variant=variant)
problem = build_riesz_map(V, grad)
for sp, expected_it in zip(solvers, expected):
> assert solve_riesz_map(problem, sp) <= expected_it
tests/regression/test_fdm.py:147:
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
tests/regression/test_fdm.py:114: in solve_riesz_map
solver.solve()
petsc4py/PETSc/Log.pyx:188: in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
???
petsc4py/PETSc/Log.pyx:189: in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
???
firedrake/adjoint_utils/variational_solver.py:89: in wrapper
out = solve(self, **kwargs)
firedrake/variational_solver.py:324: in solve
self.snes.solve(None, work)
petsc4py/PETSc/SNES.pyx:1601: in petsc4py.PETSc.SNES.solve
???
petsc4py/PETSc/PETSc.pyx:323: in petsc4py.PETSc.PetscPythonErrorHandler
???
petsc4py/PETSc/PETSc.pyx:323: in petsc4py.PETSc.PetscPythonErrorHandler
???
petsc4py/PETSc/PETSc.pyx:323: in petsc4py.PETSc.PetscPythonErrorHandler
???
petsc4py/PETSc/PETSc.pyx:323: in petsc4py.PETSc.PetscPythonErrorHandler
???
petsc4py/PETSc/PETSc.pyx:323: in petsc4py.PETSc.PetscPythonErrorHandler
???
petsc4py/PETSc/PETSc.pyx:323: in petsc4py.PETSc.PetscPythonErrorHandler
???
petsc4py/PETSc/libpetsc4py.pyx:1324: in petsc4py.PETSc.PCSetUp_Python
???
firedrake/preconditioners/base.py:129: in setUp
super().setUp(pc)
firedrake/preconditioners/base.py:45: in setUp
self.initialize(pc)
firedrake/preconditioners/pmg.py:157: in initialize
ppc.setUp()
petsc4py/PETSc/PC.pyx:563: in petsc4py.PETSc.PC.setUp
???
petsc4py/PETSc/PETSc.pyx:323: in petsc4py.PETSc.PetscPythonErrorHandler
???
petsc4py/PETSc/PETSc.pyx:323: in petsc4py.PETSc.PetscPythonErrorHandler
???
petsc4py/PETSc/PETSc.pyx:323: in petsc4py.PETSc.PetscPythonErrorHandler
???
petsc4py/PETSc/PETSc.pyx:323: in petsc4py.PETSc.PetscPythonErrorHandler
???
petsc4py/PETSc/PETSc.pyx:323: in petsc4py.PETSc.PetscPythonErrorHandler
???
petsc4py/PETSc/PETSc.pyx:323: in petsc4py.PETSc.PetscPythonErrorHandler
???
petsc4py/PETSc/PETSc.pyx:323: in petsc4py.PETSc.PetscPythonErrorHandler
???
petsc4py/PETSc/PETSc.pyx:323: in petsc4py.PETSc.PetscPythonErrorHandler
???
petsc4py/PETSc/libpetsc4py.pyx:1324: in petsc4py.PETSc.PCSetUp_Python
???
firedrake/preconditioners/base.py:129: in setUp
super().setUp(pc)
firedrake/preconditioners/base.py:45: in setUp
self.initialize(pc)
petsc4py/PETSc/Log.pyx:188: in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
???
petsc4py/PETSc/Log.pyx:189: in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
???
firedrake/preconditioners/fdm.py:178: in initialize
Amat, Pmat, self.assembly_callables = self.allocate_matrix(Amat, V_fdm, J_fdm, bcs_fdm, fcp,
petsc4py/PETSc/Log.pyx:188: in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
???
petsc4py/PETSc/Log.pyx:189: in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
???
firedrake/preconditioners/fdm.py:235: in allocate_matrix
self.lgmaps = {Vsub: Vsub.local_to_global_map([bc for bc in bcs if bc.function_space() == Vsub]) for Vsub in V}
petsc4py/PETSc/Log.pyx:188: in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
???
petsc4py/PETSc/Log.pyx:189: in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
???
firedrake/functionspaceimpl.py:840: in local_to_global_map
nodes.append(bc.nodes)
../firedrake_venv/src/PyOP2/pyop2/utils.py:62: in __get__
obj.__dict__[self.__name__] = result = self.fget(obj)
firedrake/bcs.py:169: in nodes
bcnodes1.append(hermite_stride(self._function_space.boundary_nodes(ss)))
firedrake/bcs.py:143: in hermite_stride
deriv_nodes = [i for i, node in enumerate(fe.fiat_equivalent.dual) if len(node.deriv_dict) != 0]
../firedrake_venv/src/tsfc/gem/utils.py:19: in __get__
obj.__dict__[self.__name__] = result = self.fget(obj)
../firedrake_venv/src/FInAT/finat/tensor_product.py:82: in fiat_equivalent
return FIAT.TensorProductElement(A.fiat_equivalent, B.fiat_equivalent)
../firedrake_venv/src/tsfc/gem/utils.py:19: in __get__
obj.__dict__[self.__name__] = result = self.fget(obj)
../firedrake_venv/src/FInAT/finat/cube.py:65: in fiat_equivalent
return FIAT_FlattenedDimensions(self.product.fiat_equivalent)
../firedrake_venv/src/tsfc/gem/utils.py:19: in __get__
obj.__dict__[self.__name__] = result = self.fget(obj)
../firedrake_venv/src/FInAT/finat/tensor_product.py:82: in fiat_equivalent
return FIAT.TensorProductElement(A.fiat_equivalent, B.fiat_equivalent)
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
self = <FIAT.tensor_product.TensorProductElement object at 0x7fdfb235ac60>
A = <FIAT.fdm_element.FDMLagrange object at 0x7fdfb3469160>
B = <FIAT.fdm_element.FDMLagrange object at 0x7fdfb3469160>
def __init__(self, A, B):
# set up simple things
order = min(A.get_order(), B.get_order())
if A.get_formdegree() is None or B.get_formdegree() is None:
formdegree = None
else:
formdegree = A.get_formdegree() + B.get_formdegree()
# set up reference element
ref_el = TensorProductCell(A.get_reference_element(),
B.get_reference_element())
if A.mapping()[0] != "affine" and B.mapping()[0] == "affine":
mapping = A.mapping()[0]
elif B.mapping()[0] != "affine" and A.mapping()[0] == "affine":
mapping = B.mapping()[0]
elif A.mapping()[0] == "affine" and B.mapping()[0] == "affine":
mapping = "affine"
else:
raise ValueError("check tensor product mappings - at least one must be affine")
# set up entity_ids
Adofs = A.entity_dofs()
Bdofs = B.entity_dofs()
Bsdim = B.space_dimension()
entity_ids = {}
for curAdim in Adofs:
for curBdim in Bdofs:
entity_ids[(curAdim, curBdim)] = {}
dim_cur = 0
for entityA in Adofs[curAdim]:
for entityB in Bdofs[curBdim]:
entity_ids[(curAdim, curBdim)][dim_cur] = \
[x*Bsdim + y for x in Adofs[curAdim][entityA]
for y in Bdofs[curBdim][entityB]]
dim_cur += 1
# set up dual basis
Anodes = A.dual_basis()
Bnodes = B.dual_basis()
# build the dual set by inspecting the current dual
# sets item by item.
# Currently supported cases:
# PointEval x PointEval = PointEval [scalar x scalar = scalar]
# PointScaledNormalEval x PointEval = PointScaledNormalEval [vector x scalar = vector]
# ComponentPointEvaluation x PointEval [vector x scalar = vector]
nodes = []
for Anode in Anodes:
if isinstance(Anode, functional.PointEvaluation):
for Bnode in Bnodes:
if isinstance(Bnode, functional.PointEvaluation):
# case: PointEval x PointEval
# the PointEval functional just requires the
# coordinates. these are currently stored as
# the key of a one-item dictionary. we retrieve
# these by calling get_point_dict(), and
# use the concatenation to make a new PointEval
nodes.append(functional.PointEvaluation(ref_el, _first_point(Anode) + _first_point(Bnode)))
elif isinstance(Bnode, functional.IntegralMoment):
# dummy functional for product with integral moments
nodes.append(functional.Functional(None, None, None,
{}, "Undefined"))
elif isinstance(Bnode, functional.PointDerivative):
# dummy functional for product with point derivative
nodes.append(functional.Functional(None, None, None,
{}, "Undefined"))
else:
raise NotImplementedError("unsupported functional type")
elif isinstance(Anode, functional.PointScaledNormalEvaluation):
for Bnode in Bnodes:
if isinstance(Bnode, functional.PointEvaluation):
# case: PointScaledNormalEval x PointEval
# this could be wrong if the second shape
# has spatial dimension >1, since we are not
# explicitly scaling by facet size
if len(_first_point(Bnode)) > 1:
# TODO: support this case one day
raise NotImplementedError("PointScaledNormalEval x PointEval is not yet supported if the second shape has dimension > 1")
# We cannot make a new functional.PSNEval in
# the natural way, since it tries to compute
# the normal vector by itself.
# Instead, we create things manually, and
# call Functional() with these arguments
sd = ref_el.get_spatial_dimension()
# The pt_dict is a one-item dictionary containing
# the details of the functional.
# The key is the spatial coordinate, which
# is just a concatenation of the two parts.
# The value is a list of tuples, representing
# the normal vector (scaled by the volume of
# the facet) at that point.
# Each tuple looks like (foo, (i,)); the i'th
# component of the scaled normal is foo.
# The following line is only valid when the second
# shape has spatial dimension 1 (enforced above)
Apoint, Avalue = _first_point_pair(Anode)
pt_dict = {Apoint + _first_point(Bnode): Avalue + [(0.0, (len(Apoint),))]}
# The following line should be used in the
# general case
# pt_dict = {Anode.get_point_dict().keys()[0] + Bnode.get_point_dict().keys()[0]: Anode.get_point_dict().values()[0] + [(0.0, (ii,)) for ii in range(len(Anode.get_point_dict().keys()[0]), len(Anode.get_point_dict().keys()[0]) + len(Bnode.get_point_dict().keys()[0]))]}
# THE FOLLOWING IS PROBABLY CORRECT BUT UNTESTED
shp = (sd,)
nodes.append(functional.Functional(ref_el, shp, pt_dict, {}, "PointScaledNormalEval"))
else:
raise NotImplementedError("unsupported functional type")
elif isinstance(Anode, functional.PointEdgeTangentEvaluation):
for Bnode in Bnodes:
if isinstance(Bnode, functional.PointEvaluation):
# case: PointEdgeTangentEval x PointEval
# this is very similar to the case above, so comments omitted
if len(_first_point(Bnode)) > 1:
raise NotImplementedError("PointEdgeTangentEval x PointEval is not yet supported if the second shape has dimension > 1")
sd = ref_el.get_spatial_dimension()
Apoint, Avalue = _first_point_pair(Anode)
pt_dict = {Apoint + _first_point(Bnode): Avalue + [(0.0, (len(Apoint),))]}
# THE FOLLOWING IS PROBABLY CORRECT BUT UNTESTED
shp = (sd,)
nodes.append(functional.Functional(ref_el, shp, pt_dict, {}, "PointEdgeTangent"))
else:
raise NotImplementedError("unsupported functional type")
elif isinstance(Anode, functional.ComponentPointEvaluation):
for Bnode in Bnodes:
if isinstance(Bnode, functional.PointEvaluation):
# case: ComponentPointEval x PointEval
# the CptPointEval functional requires the component
# and the coordinates. very similar to PE x PE case.
sd = ref_el.get_spatial_dimension()
nodes.append(functional.ComponentPointEvaluation(ref_el, Anode.comp, (sd,), _first_point(Anode) + _first_point(Bnode)))
else:
raise NotImplementedError("unsupported functional type")
elif isinstance(Anode, functional.FrobeniusIntegralMoment):
for Bnode in Bnodes:
if isinstance(Bnode, functional.PointEvaluation):
# case: FroIntMom x PointEval
sd = ref_el.get_spatial_dimension()
pt_dict = {}
pt_old = Anode.get_point_dict()
for pt in pt_old:
pt_dict[pt+_first_point(Bnode)] = pt_old[pt] + [(0.0, sd-1)]
# THE FOLLOWING IS PROBABLY CORRECT BUT UNTESTED
shp = (sd,)
nodes.append(functional.Functional(ref_el, shp, pt_dict, {}, "FrobeniusIntegralMoment"))
else:
raise NotImplementedError("unsupported functional type")
elif isinstance(Anode, functional.IntegralMoment):
for Bnode in Bnodes:
if isinstance(Bnode, functional.PointEvaluation):
# case: IntMom x PointEval
sd = ref_el.get_spatial_dimension()
pt_dict = {}
pt_old = Anode.get_point_dict()
for pt in pt_old:
pt_dict[pt+_first_point(Bnode)] = pt_old[pt]
# THE FOLLOWING IS PROBABLY CORRECT BUT UNTESTED
shp = (sd,)
nodes.append(functional.Functional(ref_el, shp, pt_dict, {}, "IntegralMoment"))
else:
> raise NotImplementedError("unsupported functional type")
E NotImplementedError: unsupported functional type
../firedrake_venv/src/fiat/FIAT/tensor_product.py:198: NotImplementedError
Check failure on line 114 in tests/regression/test_fdm.py
github-actions / Firedrake real
test_fdm.test_p_independence_hgrad[Box-fdm]
NotImplementedError: unsupported functional type
Raw output
mesh = Mesh(VectorElement(TensorProductElement(FiniteElement('Q', quadrilateral, 1), FiniteElement('Lagrange', interval, 1), cell=TensorProductCell(quadrilateral, interval)), dim=3), 3257)
variant = 'fdm'
@pytest.mark.skipcomplex
def test_p_independence_hgrad(mesh, variant):
family = "Lagrange"
expected = [16, 12] if mesh.topological_dimension() == 3 else [9, 7]
solvers = [fdmstar] if variant is None else [fdmstar, facetstar]
for degree in range(3, 6):
V = FunctionSpace(mesh, family, degree, variant=variant)
problem = build_riesz_map(V, grad)
for sp, expected_it in zip(solvers, expected):
> assert solve_riesz_map(problem, sp) <= expected_it
tests/regression/test_fdm.py:147:
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
tests/regression/test_fdm.py:114: in solve_riesz_map
solver.solve()
petsc4py/PETSc/Log.pyx:188: in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
???
petsc4py/PETSc/Log.pyx:189: in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
???
firedrake/adjoint_utils/variational_solver.py:89: in wrapper
out = solve(self, **kwargs)
firedrake/variational_solver.py:307: in solve
dbc.apply(problem.u_restrict)
petsc4py/PETSc/Log.pyx:188: in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
???
petsc4py/PETSc/Log.pyx:189: in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
???
firedrake/adjoint_utils/dirichletbc.py:32: in wrapper
ret = apply(self, *args, **kwargs)
firedrake/bcs.py:451: in apply
r.assign(self.function_arg, subset=self.node_set)
../firedrake_venv/src/PyOP2/pyop2/utils.py:62: in __get__
obj.__dict__[self.__name__] = result = self.fget(obj)
firedrake/bcs.py:179: in node_set
return op2.Subset(self._function_space.node_set, self.nodes)
../firedrake_venv/src/PyOP2/pyop2/utils.py:62: in __get__
obj.__dict__[self.__name__] = result = self.fget(obj)
firedrake/bcs.py:169: in nodes
bcnodes1.append(hermite_stride(self._function_space.boundary_nodes(ss)))
firedrake/bcs.py:143: in hermite_stride
deriv_nodes = [i for i, node in enumerate(fe.fiat_equivalent.dual) if len(node.deriv_dict) != 0]
../firedrake_venv/src/tsfc/gem/utils.py:19: in __get__
obj.__dict__[self.__name__] = result = self.fget(obj)
../firedrake_venv/src/FInAT/finat/tensor_product.py:82: in fiat_equivalent
return FIAT.TensorProductElement(A.fiat_equivalent, B.fiat_equivalent)
../firedrake_venv/src/tsfc/gem/utils.py:19: in __get__
obj.__dict__[self.__name__] = result = self.fget(obj)
../firedrake_venv/src/FInAT/finat/cube.py:65: in fiat_equivalent
return FIAT_FlattenedDimensions(self.product.fiat_equivalent)
../firedrake_venv/src/tsfc/gem/utils.py:19: in __get__
obj.__dict__[self.__name__] = result = self.fget(obj)
../firedrake_venv/src/FInAT/finat/tensor_product.py:82: in fiat_equivalent
return FIAT.TensorProductElement(A.fiat_equivalent, B.fiat_equivalent)
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
self = <FIAT.tensor_product.TensorProductElement object at 0x7fdfb23a3170>
A = <FIAT.fdm_element.FDMLagrange object at 0x7fdfb3469160>
B = <FIAT.fdm_element.FDMLagrange object at 0x7fdfb3469160>
def __init__(self, A, B):
# set up simple things
order = min(A.get_order(), B.get_order())
if A.get_formdegree() is None or B.get_formdegree() is None:
formdegree = None
else:
formdegree = A.get_formdegree() + B.get_formdegree()
# set up reference element
ref_el = TensorProductCell(A.get_reference_element(),
B.get_reference_element())
if A.mapping()[0] != "affine" and B.mapping()[0] == "affine":
mapping = A.mapping()[0]
elif B.mapping()[0] != "affine" and A.mapping()[0] == "affine":
mapping = B.mapping()[0]
elif A.mapping()[0] == "affine" and B.mapping()[0] == "affine":
mapping = "affine"
else:
raise ValueError("check tensor product mappings - at least one must be affine")
# set up entity_ids
Adofs = A.entity_dofs()
Bdofs = B.entity_dofs()
Bsdim = B.space_dimension()
entity_ids = {}
for curAdim in Adofs:
for curBdim in Bdofs:
entity_ids[(curAdim, curBdim)] = {}
dim_cur = 0
for entityA in Adofs[curAdim]:
for entityB in Bdofs[curBdim]:
entity_ids[(curAdim, curBdim)][dim_cur] = \
[x*Bsdim + y for x in Adofs[curAdim][entityA]
for y in Bdofs[curBdim][entityB]]
dim_cur += 1
# set up dual basis
Anodes = A.dual_basis()
Bnodes = B.dual_basis()
# build the dual set by inspecting the current dual
# sets item by item.
# Currently supported cases:
# PointEval x PointEval = PointEval [scalar x scalar = scalar]
# PointScaledNormalEval x PointEval = PointScaledNormalEval [vector x scalar = vector]
# ComponentPointEvaluation x PointEval [vector x scalar = vector]
nodes = []
for Anode in Anodes:
if isinstance(Anode, functional.PointEvaluation):
for Bnode in Bnodes:
if isinstance(Bnode, functional.PointEvaluation):
# case: PointEval x PointEval
# the PointEval functional just requires the
# coordinates. these are currently stored as
# the key of a one-item dictionary. we retrieve
# these by calling get_point_dict(), and
# use the concatenation to make a new PointEval
nodes.append(functional.PointEvaluation(ref_el, _first_point(Anode) + _first_point(Bnode)))
elif isinstance(Bnode, functional.IntegralMoment):
# dummy functional for product with integral moments
nodes.append(functional.Functional(None, None, None,
{}, "Undefined"))
elif isinstance(Bnode, functional.PointDerivative):
# dummy functional for product with point derivative
nodes.append(functional.Functional(None, None, None,
{}, "Undefined"))
else:
raise NotImplementedError("unsupported functional type")
elif isinstance(Anode, functional.PointScaledNormalEvaluation):
for Bnode in Bnodes:
if isinstance(Bnode, functional.PointEvaluation):
# case: PointScaledNormalEval x PointEval
# this could be wrong if the second shape
# has spatial dimension >1, since we are not
# explicitly scaling by facet size
if len(_first_point(Bnode)) > 1:
# TODO: support this case one day
raise NotImplementedError("PointScaledNormalEval x PointEval is not yet supported if the second shape has dimension > 1")
# We cannot make a new functional.PSNEval in
# the natural way, since it tries to compute
# the normal vector by itself.
# Instead, we create things manually, and
# call Functional() with these arguments
sd = ref_el.get_spatial_dimension()
# The pt_dict is a one-item dictionary containing
# the details of the functional.
# The key is the spatial coordinate, which
# is just a concatenation of the two parts.
# The value is a list of tuples, representing
# the normal vector (scaled by the volume of
# the facet) at that point.
# Each tuple looks like (foo, (i,)); the i'th
# component of the scaled normal is foo.
# The following line is only valid when the second
# shape has spatial dimension 1 (enforced above)
Apoint, Avalue = _first_point_pair(Anode)
pt_dict = {Apoint + _first_point(Bnode): Avalue + [(0.0, (len(Apoint),))]}
# The following line should be used in the
# general case
# pt_dict = {Anode.get_point_dict().keys()[0] + Bnode.get_point_dict().keys()[0]: Anode.get_point_dict().values()[0] + [(0.0, (ii,)) for ii in range(len(Anode.get_point_dict().keys()[0]), len(Anode.get_point_dict().keys()[0]) + len(Bnode.get_point_dict().keys()[0]))]}
# THE FOLLOWING IS PROBABLY CORRECT BUT UNTESTED
shp = (sd,)
nodes.append(functional.Functional(ref_el, shp, pt_dict, {}, "PointScaledNormalEval"))
else:
raise NotImplementedError("unsupported functional type")
elif isinstance(Anode, functional.PointEdgeTangentEvaluation):
for Bnode in Bnodes:
if isinstance(Bnode, functional.PointEvaluation):
# case: PointEdgeTangentEval x PointEval
# this is very similar to the case above, so comments omitted
if len(_first_point(Bnode)) > 1:
raise NotImplementedError("PointEdgeTangentEval x PointEval is not yet supported if the second shape has dimension > 1")
sd = ref_el.get_spatial_dimension()
Apoint, Avalue = _first_point_pair(Anode)
pt_dict = {Apoint + _first_point(Bnode): Avalue + [(0.0, (len(Apoint),))]}
# THE FOLLOWING IS PROBABLY CORRECT BUT UNTESTED
shp = (sd,)
nodes.append(functional.Functional(ref_el, shp, pt_dict, {}, "PointEdgeTangent"))
else:
raise NotImplementedError("unsupported functional type")
elif isinstance(Anode, functional.ComponentPointEvaluation):
for Bnode in Bnodes:
if isinstance(Bnode, functional.PointEvaluation):
# case: ComponentPointEval x PointEval
# the CptPointEval functional requires the component
# and the coordinates. very similar to PE x PE case.
sd = ref_el.get_spatial_dimension()
nodes.append(functional.ComponentPointEvaluation(ref_el, Anode.comp, (sd,), _first_point(Anode) + _first_point(Bnode)))
else:
raise NotImplementedError("unsupported functional type")
elif isinstance(Anode, functional.FrobeniusIntegralMoment):
for Bnode in Bnodes:
if isinstance(Bnode, functional.PointEvaluation):
# case: FroIntMom x PointEval
sd = ref_el.get_spatial_dimension()
pt_dict = {}
pt_old = Anode.get_point_dict()
for pt in pt_old:
pt_dict[pt+_first_point(Bnode)] = pt_old[pt] + [(0.0, sd-1)]
# THE FOLLOWING IS PROBABLY CORRECT BUT UNTESTED
shp = (sd,)
nodes.append(functional.Functional(ref_el, shp, pt_dict, {}, "FrobeniusIntegralMoment"))
else:
raise NotImplementedError("unsupported functional type")
elif isinstance(Anode, functional.IntegralMoment):
for Bnode in Bnodes:
if isinstance(Bnode, functional.PointEvaluation):
# case: IntMom x PointEval
sd = ref_el.get_spatial_dimension()
pt_dict = {}
pt_old = Anode.get_point_dict()
for pt in pt_old:
pt_dict[pt+_first_point(Bnode)] = pt_old[pt]
# THE FOLLOWING IS PROBABLY CORRECT BUT UNTESTED
shp = (sd,)
nodes.append(functional.Functional(ref_el, shp, pt_dict, {}, "IntegralMoment"))
else:
> raise NotImplementedError("unsupported functional type")
E NotImplementedError: unsupported functional type
../firedrake_venv/src/fiat/FIAT/tensor_product.py:198: NotImplementedError
Check failure on line 203 in tests/regression/test_fdm.py
github-actions / Firedrake real
test_fdm.test_variable_coefficient[Rectangle]
NotImplementedError: unsupported functional type
Raw output
mesh = Mesh(VectorElement(FiniteElement('Q', quadrilateral, 1), dim=2), 3930)
@pytest.mark.skipcomplex
def test_variable_coefficient(mesh):
gdim = mesh.geometric_dimension()
k = 4
V = FunctionSpace(mesh, "Lagrange", k)
u = TrialFunction(V)
v = TestFunction(V)
x = SpatialCoordinate(mesh)
x -= Constant([0.5]*len(x))
# variable coefficients
alphas = [0.1+10*dot(x, x)]*gdim
alphas[0] = 1+10*exp(-dot(x, x))
alpha = diag(as_vector(alphas))
beta = ((10*cos(3*pi*x[0]) + 20*sin(2*pi*x[1]))*cos(pi*x[gdim-1]))**2
a = (inner(grad(v), dot(alpha, grad(u))) + inner(v, beta*u))*dx(degree=3*k+2)
L = inner(v, Constant(1))*dx
subs = ("on_boundary",)
if mesh.cell_set._extruded:
subs += ("top", "bottom")
bcs = [DirichletBC(V, zero(V.ufl_element().value_shape), sub) for sub in subs]
uh = Function(V)
problem = LinearVariationalProblem(a, L, uh, bcs=bcs)
solver = LinearVariationalSolver(problem, solver_parameters=fdmstar)
> solver.solve()
tests/regression/test_fdm.py:203:
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
petsc4py/PETSc/Log.pyx:188: in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
???
petsc4py/PETSc/Log.pyx:189: in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
???
firedrake/adjoint_utils/variational_solver.py:89: in wrapper
out = solve(self, **kwargs)
firedrake/variational_solver.py:324: in solve
self.snes.solve(None, work)
petsc4py/PETSc/SNES.pyx:1601: in petsc4py.PETSc.SNES.solve
???
petsc4py/PETSc/PETSc.pyx:323: in petsc4py.PETSc.PetscPythonErrorHandler
???
petsc4py/PETSc/PETSc.pyx:323: in petsc4py.PETSc.PetscPythonErrorHandler
???
petsc4py/PETSc/PETSc.pyx:323: in petsc4py.PETSc.PetscPythonErrorHandler
???
petsc4py/PETSc/PETSc.pyx:323: in petsc4py.PETSc.PetscPythonErrorHandler
???
petsc4py/PETSc/PETSc.pyx:323: in petsc4py.PETSc.PetscPythonErrorHandler
???
petsc4py/PETSc/PETSc.pyx:323: in petsc4py.PETSc.PetscPythonErrorHandler
???
petsc4py/PETSc/libpetsc4py.pyx:1324: in petsc4py.PETSc.PCSetUp_Python
???
firedrake/preconditioners/base.py:129: in setUp
super().setUp(pc)
firedrake/preconditioners/base.py:45: in setUp
self.initialize(pc)
firedrake/preconditioners/pmg.py:157: in initialize
ppc.setUp()
petsc4py/PETSc/PC.pyx:563: in petsc4py.PETSc.PC.setUp
???
petsc4py/PETSc/PETSc.pyx:323: in petsc4py.PETSc.PetscPythonErrorHandler
???
petsc4py/PETSc/PETSc.pyx:323: in petsc4py.PETSc.PetscPythonErrorHandler
???
petsc4py/PETSc/PETSc.pyx:323: in petsc4py.PETSc.PetscPythonErrorHandler
???
petsc4py/PETSc/PETSc.pyx:323: in petsc4py.PETSc.PetscPythonErrorHandler
???
petsc4py/PETSc/PETSc.pyx:323: in petsc4py.PETSc.PetscPythonErrorHandler
???
petsc4py/PETSc/PETSc.pyx:323: in petsc4py.PETSc.PetscPythonErrorHandler
???
petsc4py/PETSc/PETSc.pyx:323: in petsc4py.PETSc.PetscPythonErrorHandler
???
petsc4py/PETSc/PETSc.pyx:323: in petsc4py.PETSc.PetscPythonErrorHandler
???
petsc4py/PETSc/libpetsc4py.pyx:1324: in petsc4py.PETSc.PCSetUp_Python
???
firedrake/preconditioners/base.py:129: in setUp
super().setUp(pc)
firedrake/preconditioners/base.py:45: in setUp
self.initialize(pc)
petsc4py/PETSc/Log.pyx:188: in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
???
petsc4py/PETSc/Log.pyx:189: in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
???
firedrake/preconditioners/fdm.py:178: in initialize
Amat, Pmat, self.assembly_callables = self.allocate_matrix(Amat, V_fdm, J_fdm, bcs_fdm, fcp,
petsc4py/PETSc/Log.pyx:188: in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
???
petsc4py/PETSc/Log.pyx:189: in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
???
firedrake/preconditioners/fdm.py:235: in allocate_matrix
self.lgmaps = {Vsub: Vsub.local_to_global_map([bc for bc in bcs if bc.function_space() == Vsub]) for Vsub in V}
petsc4py/PETSc/Log.pyx:188: in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
???
petsc4py/PETSc/Log.pyx:189: in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
???
firedrake/functionspaceimpl.py:840: in local_to_global_map
nodes.append(bc.nodes)
../firedrake_venv/src/PyOP2/pyop2/utils.py:62: in __get__
obj.__dict__[self.__name__] = result = self.fget(obj)
firedrake/bcs.py:154: in nodes
bcnodes.append(hermite_stride(self._function_space.boundary_nodes(s)))
firedrake/bcs.py:143: in hermite_stride
deriv_nodes = [i for i, node in enumerate(fe.fiat_equivalent.dual) if len(node.deriv_dict) != 0]
../firedrake_venv/src/tsfc/gem/utils.py:19: in __get__
obj.__dict__[self.__name__] = result = self.fget(obj)
../firedrake_venv/src/FInAT/finat/cube.py:65: in fiat_equivalent
return FIAT_FlattenedDimensions(self.product.fiat_equivalent)
../firedrake_venv/src/tsfc/gem/utils.py:19: in __get__
obj.__dict__[self.__name__] = result = self.fget(obj)
../firedrake_venv/src/FInAT/finat/tensor_product.py:82: in fiat_equivalent
return FIAT.TensorProductElement(A.fiat_equivalent, B.fiat_equivalent)
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
self = <FIAT.tensor_product.TensorProductElement object at 0x7fdfa82d4350>
A = <FIAT.fdm_element.FDMLagrange object at 0x7fdfb075ffb0>
B = <FIAT.fdm_element.FDMLagrange object at 0x7fdfb075ffb0>
def __init__(self, A, B):
# set up simple things
order = min(A.get_order(), B.get_order())
if A.get_formdegree() is None or B.get_formdegree() is None:
formdegree = None
else:
formdegree = A.get_formdegree() + B.get_formdegree()
# set up reference element
ref_el = TensorProductCell(A.get_reference_element(),
B.get_reference_element())
if A.mapping()[0] != "affine" and B.mapping()[0] == "affine":
mapping = A.mapping()[0]
elif B.mapping()[0] != "affine" and A.mapping()[0] == "affine":
mapping = B.mapping()[0]
elif A.mapping()[0] == "affine" and B.mapping()[0] == "affine":
mapping = "affine"
else:
raise ValueError("check tensor product mappings - at least one must be affine")
# set up entity_ids
Adofs = A.entity_dofs()
Bdofs = B.entity_dofs()
Bsdim = B.space_dimension()
entity_ids = {}
for curAdim in Adofs:
for curBdim in Bdofs:
entity_ids[(curAdim, curBdim)] = {}
dim_cur = 0
for entityA in Adofs[curAdim]:
for entityB in Bdofs[curBdim]:
entity_ids[(curAdim, curBdim)][dim_cur] = \
[x*Bsdim + y for x in Adofs[curAdim][entityA]
for y in Bdofs[curBdim][entityB]]
dim_cur += 1
# set up dual basis
Anodes = A.dual_basis()
Bnodes = B.dual_basis()
# build the dual set by inspecting the current dual
# sets item by item.
# Currently supported cases:
# PointEval x PointEval = PointEval [scalar x scalar = scalar]
# PointScaledNormalEval x PointEval = PointScaledNormalEval [vector x scalar = vector]
# ComponentPointEvaluation x PointEval [vector x scalar = vector]
nodes = []
for Anode in Anodes:
if isinstance(Anode, functional.PointEvaluation):
for Bnode in Bnodes:
if isinstance(Bnode, functional.PointEvaluation):
# case: PointEval x PointEval
# the PointEval functional just requires the
# coordinates. these are currently stored as
# the key of a one-item dictionary. we retrieve
# these by calling get_point_dict(), and
# use the concatenation to make a new PointEval
nodes.append(functional.PointEvaluation(ref_el, _first_point(Anode) + _first_point(Bnode)))
elif isinstance(Bnode, functional.IntegralMoment):
# dummy functional for product with integral moments
nodes.append(functional.Functional(None, None, None,
{}, "Undefined"))
elif isinstance(Bnode, functional.PointDerivative):
# dummy functional for product with point derivative
nodes.append(functional.Functional(None, None, None,
{}, "Undefined"))
else:
raise NotImplementedError("unsupported functional type")
elif isinstance(Anode, functional.PointScaledNormalEvaluation):
for Bnode in Bnodes:
if isinstance(Bnode, functional.PointEvaluation):
# case: PointScaledNormalEval x PointEval
# this could be wrong if the second shape
# has spatial dimension >1, since we are not
# explicitly scaling by facet size
if len(_first_point(Bnode)) > 1:
# TODO: support this case one day
raise NotImplementedError("PointScaledNormalEval x PointEval is not yet supported if the second shape has dimension > 1")
# We cannot make a new functional.PSNEval in
# the natural way, since it tries to compute
# the normal vector by itself.
# Instead, we create things manually, and
# call Functional() with these arguments
sd = ref_el.get_spatial_dimension()
# The pt_dict is a one-item dictionary containing
# the details of the functional.
# The key is the spatial coordinate, which
# is just a concatenation of the two parts.
# The value is a list of tuples, representing
# the normal vector (scaled by the volume of
# the facet) at that point.
# Each tuple looks like (foo, (i,)); the i'th
# component of the scaled normal is foo.
# The following line is only valid when the second
# shape has spatial dimension 1 (enforced above)
Apoint, Avalue = _first_point_pair(Anode)
pt_dict = {Apoint + _first_point(Bnode): Avalue + [(0.0, (len(Apoint),))]}
# The following line should be used in the
# general case
# pt_dict = {Anode.get_point_dict().keys()[0] + Bnode.get_point_dict().keys()[0]: Anode.get_point_dict().values()[0] + [(0.0, (ii,)) for ii in range(len(Anode.get_point_dict().keys()[0]), len(Anode.get_point_dict().keys()[0]) + len(Bnode.get_point_dict().keys()[0]))]}
# THE FOLLOWING IS PROBABLY CORRECT BUT UNTESTED
shp = (sd,)
nodes.append(functional.Functional(ref_el, shp, pt_dict, {}, "PointScaledNormalEval"))
else:
raise NotImplementedError("unsupported functional type")
elif isinstance(Anode, functional.PointEdgeTangentEvaluation):
for Bnode in Bnodes:
if isinstance(Bnode, functional.PointEvaluation):
# case: PointEdgeTangentEval x PointEval
# this is very similar to the case above, so comments omitted
if len(_first_point(Bnode)) > 1:
raise NotImplementedError("PointEdgeTangentEval x PointEval is not yet supported if the second shape has dimension > 1")
sd = ref_el.get_spatial_dimension()
Apoint, Avalue = _first_point_pair(Anode)
pt_dict = {Apoint + _first_point(Bnode): Avalue + [(0.0, (len(Apoint),))]}
# THE FOLLOWING IS PROBABLY CORRECT BUT UNTESTED
shp = (sd,)
nodes.append(functional.Functional(ref_el, shp, pt_dict, {}, "PointEdgeTangent"))
else:
raise NotImplementedError("unsupported functional type")
elif isinstance(Anode, functional.ComponentPointEvaluation):
for Bnode in Bnodes:
if isinstance(Bnode, functional.PointEvaluation):
# case: ComponentPointEval x PointEval
# the CptPointEval functional requires the component
# and the coordinates. very similar to PE x PE case.
sd = ref_el.get_spatial_dimension()
nodes.append(functional.ComponentPointEvaluation(ref_el, Anode.comp, (sd,), _first_point(Anode) + _first_point(Bnode)))
else:
raise NotImplementedError("unsupported functional type")
elif isinstance(Anode, functional.FrobeniusIntegralMoment):
for Bnode in Bnodes:
if isinstance(Bnode, functional.PointEvaluation):
# case: FroIntMom x PointEval
sd = ref_el.get_spatial_dimension()
pt_dict = {}
pt_old = Anode.get_point_dict()
for pt in pt_old:
pt_dict[pt+_first_point(Bnode)] = pt_old[pt] + [(0.0, sd-1)]
# THE FOLLOWING IS PROBABLY CORRECT BUT UNTESTED
shp = (sd,)
nodes.append(functional.Functional(ref_el, shp, pt_dict, {}, "FrobeniusIntegralMoment"))
else:
raise NotImplementedError("unsupported functional type")
elif isinstance(Anode, functional.IntegralMoment):
for Bnode in Bnodes:
if isinstance(Bnode, functional.PointEvaluation):
# case: IntMom x PointEval
sd = ref_el.get_spatial_dimension()
pt_dict = {}
pt_old = Anode.get_point_dict()
for pt in pt_old:
pt_dict[pt+_first_point(Bnode)] = pt_old[pt]
# THE FOLLOWING IS PROBABLY CORRECT BUT UNTESTED
shp = (sd,)
nodes.append(functional.Functional(ref_el, shp, pt_dict, {}, "IntegralMoment"))
else:
> raise NotImplementedError("unsupported functional type")
E NotImplementedError: unsupported functional type
../firedrake_venv/src/fiat/FIAT/tensor_product.py:198: NotImplementedError
Check failure on line 203 in tests/regression/test_fdm.py
github-actions / Firedrake real
test_fdm.test_variable_coefficient[Box]
NotImplementedError: unsupported functional type
Raw output
mesh = Mesh(VectorElement(TensorProductElement(FiniteElement('Q', quadrilateral, 1), FiniteElement('Lagrange', interval, 1), cell=TensorProductCell(quadrilateral, interval)), dim=3), 3958)
@pytest.mark.skipcomplex
def test_variable_coefficient(mesh):
gdim = mesh.geometric_dimension()
k = 4
V = FunctionSpace(mesh, "Lagrange", k)
u = TrialFunction(V)
v = TestFunction(V)
x = SpatialCoordinate(mesh)
x -= Constant([0.5]*len(x))
# variable coefficients
alphas = [0.1+10*dot(x, x)]*gdim
alphas[0] = 1+10*exp(-dot(x, x))
alpha = diag(as_vector(alphas))
beta = ((10*cos(3*pi*x[0]) + 20*sin(2*pi*x[1]))*cos(pi*x[gdim-1]))**2
a = (inner(grad(v), dot(alpha, grad(u))) + inner(v, beta*u))*dx(degree=3*k+2)
L = inner(v, Constant(1))*dx
subs = ("on_boundary",)
if mesh.cell_set._extruded:
subs += ("top", "bottom")
bcs = [DirichletBC(V, zero(V.ufl_element().value_shape), sub) for sub in subs]
uh = Function(V)
problem = LinearVariationalProblem(a, L, uh, bcs=bcs)
solver = LinearVariationalSolver(problem, solver_parameters=fdmstar)
> solver.solve()
tests/regression/test_fdm.py:203:
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
petsc4py/PETSc/Log.pyx:188: in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
???
petsc4py/PETSc/Log.pyx:189: in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
???
firedrake/adjoint_utils/variational_solver.py:89: in wrapper
out = solve(self, **kwargs)
firedrake/variational_solver.py:324: in solve
self.snes.solve(None, work)
petsc4py/PETSc/SNES.pyx:1601: in petsc4py.PETSc.SNES.solve
???
petsc4py/PETSc/PETSc.pyx:323: in petsc4py.PETSc.PetscPythonErrorHandler
???
petsc4py/PETSc/PETSc.pyx:323: in petsc4py.PETSc.PetscPythonErrorHandler
???
petsc4py/PETSc/PETSc.pyx:323: in petsc4py.PETSc.PetscPythonErrorHandler
???
petsc4py/PETSc/PETSc.pyx:323: in petsc4py.PETSc.PetscPythonErrorHandler
???
petsc4py/PETSc/PETSc.pyx:323: in petsc4py.PETSc.PetscPythonErrorHandler
???
petsc4py/PETSc/PETSc.pyx:323: in petsc4py.PETSc.PetscPythonErrorHandler
???
petsc4py/PETSc/libpetsc4py.pyx:1324: in petsc4py.PETSc.PCSetUp_Python
???
firedrake/preconditioners/base.py:129: in setUp
super().setUp(pc)
firedrake/preconditioners/base.py:45: in setUp
self.initialize(pc)
firedrake/preconditioners/pmg.py:157: in initialize
ppc.setUp()
petsc4py/PETSc/PC.pyx:563: in petsc4py.PETSc.PC.setUp
???
petsc4py/PETSc/PETSc.pyx:323: in petsc4py.PETSc.PetscPythonErrorHandler
???
petsc4py/PETSc/PETSc.pyx:323: in petsc4py.PETSc.PetscPythonErrorHandler
???
petsc4py/PETSc/PETSc.pyx:323: in petsc4py.PETSc.PetscPythonErrorHandler
???
petsc4py/PETSc/PETSc.pyx:323: in petsc4py.PETSc.PetscPythonErrorHandler
???
petsc4py/PETSc/PETSc.pyx:323: in petsc4py.PETSc.PetscPythonErrorHandler
???
petsc4py/PETSc/PETSc.pyx:323: in petsc4py.PETSc.PetscPythonErrorHandler
???
petsc4py/PETSc/PETSc.pyx:323: in petsc4py.PETSc.PetscPythonErrorHandler
???
petsc4py/PETSc/PETSc.pyx:323: in petsc4py.PETSc.PetscPythonErrorHandler
???
petsc4py/PETSc/libpetsc4py.pyx:1324: in petsc4py.PETSc.PCSetUp_Python
???
firedrake/preconditioners/base.py:129: in setUp
super().setUp(pc)
firedrake/preconditioners/base.py:45: in setUp
self.initialize(pc)
petsc4py/PETSc/Log.pyx:188: in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
???
petsc4py/PETSc/Log.pyx:189: in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
???
firedrake/preconditioners/fdm.py:178: in initialize
Amat, Pmat, self.assembly_callables = self.allocate_matrix(Amat, V_fdm, J_fdm, bcs_fdm, fcp,
petsc4py/PETSc/Log.pyx:188: in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
???
petsc4py/PETSc/Log.pyx:189: in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
???
firedrake/preconditioners/fdm.py:235: in allocate_matrix
self.lgmaps = {Vsub: Vsub.local_to_global_map([bc for bc in bcs if bc.function_space() == Vsub]) for Vsub in V}
petsc4py/PETSc/Log.pyx:188: in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
???
petsc4py/PETSc/Log.pyx:189: in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
???
firedrake/functionspaceimpl.py:840: in local_to_global_map
nodes.append(bc.nodes)
../firedrake_venv/src/PyOP2/pyop2/utils.py:62: in __get__
obj.__dict__[self.__name__] = result = self.fget(obj)
firedrake/bcs.py:154: in nodes
bcnodes.append(hermite_stride(self._function_space.boundary_nodes(s)))
firedrake/bcs.py:143: in hermite_stride
deriv_nodes = [i for i, node in enumerate(fe.fiat_equivalent.dual) if len(node.deriv_dict) != 0]
../firedrake_venv/src/tsfc/gem/utils.py:19: in __get__
obj.__dict__[self.__name__] = result = self.fget(obj)
../firedrake_venv/src/FInAT/finat/tensor_product.py:82: in fiat_equivalent
return FIAT.TensorProductElement(A.fiat_equivalent, B.fiat_equivalent)
../firedrake_venv/src/tsfc/gem/utils.py:19: in __get__
obj.__dict__[self.__name__] = result = self.fget(obj)
../firedrake_venv/src/FInAT/finat/cube.py:65: in fiat_equivalent
return FIAT_FlattenedDimensions(self.product.fiat_equivalent)
../firedrake_venv/src/tsfc/gem/utils.py:19: in __get__
obj.__dict__[self.__name__] = result = self.fget(obj)
../firedrake_venv/src/FInAT/finat/tensor_product.py:82: in fiat_equivalent
return FIAT.TensorProductElement(A.fiat_equivalent, B.fiat_equivalent)
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
self = <FIAT.tensor_product.TensorProductElement object at 0x7fdfa82d2f30>
A = <FIAT.fdm_element.FDMLagrange object at 0x7fdfb075ffb0>
B = <FIAT.fdm_element.FDMLagrange object at 0x7fdfb075ffb0>
def __init__(self, A, B):
# set up simple things
order = min(A.get_order(), B.get_order())
if A.get_formdegree() is None or B.get_formdegree() is None:
formdegree = None
else:
formdegree = A.get_formdegree() + B.get_formdegree()
# set up reference element
ref_el = TensorProductCell(A.get_reference_element(),
B.get_reference_element())
if A.mapping()[0] != "affine" and B.mapping()[0] == "affine":
mapping = A.mapping()[0]
elif B.mapping()[0] != "affine" and A.mapping()[0] == "affine":
mapping = B.mapping()[0]
elif A.mapping()[0] == "affine" and B.mapping()[0] == "affine":
mapping = "affine"
else:
raise ValueError("check tensor product mappings - at least one must be affine")
# set up entity_ids
Adofs = A.entity_dofs()
Bdofs = B.entity_dofs()
Bsdim = B.space_dimension()
entity_ids = {}
for curAdim in Adofs:
for curBdim in Bdofs:
entity_ids[(curAdim, curBdim)] = {}
dim_cur = 0
for entityA in Adofs[curAdim]:
for entityB in Bdofs[curBdim]:
entity_ids[(curAdim, curBdim)][dim_cur] = \
[x*Bsdim + y for x in Adofs[curAdim][entityA]
for y in Bdofs[curBdim][entityB]]
dim_cur += 1
# set up dual basis
Anodes = A.dual_basis()
Bnodes = B.dual_basis()
# build the dual set by inspecting the current dual
# sets item by item.
# Currently supported cases:
# PointEval x PointEval = PointEval [scalar x scalar = scalar]
# PointScaledNormalEval x PointEval = PointScaledNormalEval [vector x scalar = vector]
# ComponentPointEvaluation x PointEval [vector x scalar = vector]
nodes = []
for Anode in Anodes:
if isinstance(Anode, functional.PointEvaluation):
for Bnode in Bnodes:
if isinstance(Bnode, functional.PointEvaluation):
# case: PointEval x PointEval
# the PointEval functional just requires the
# coordinates. these are currently stored as
# the key of a one-item dictionary. we retrieve
# these by calling get_point_dict(), and
# use the concatenation to make a new PointEval
nodes.append(functional.PointEvaluation(ref_el, _first_point(Anode) + _first_point(Bnode)))
elif isinstance(Bnode, functional.IntegralMoment):
# dummy functional for product with integral moments
nodes.append(functional.Functional(None, None, None,
{}, "Undefined"))
elif isinstance(Bnode, functional.PointDerivative):
# dummy functional for product with point derivative
nodes.append(functional.Functional(None, None, None,
{}, "Undefined"))
else:
raise NotImplementedError("unsupported functional type")
elif isinstance(Anode, functional.PointScaledNormalEvaluation):
for Bnode in Bnodes:
if isinstance(Bnode, functional.PointEvaluation):
# case: PointScaledNormalEval x PointEval
# this could be wrong if the second shape
# has spatial dimension >1, since we are not
# explicitly scaling by facet size
if len(_first_point(Bnode)) > 1:
# TODO: support this case one day
raise NotImplementedError("PointScaledNormalEval x PointEval is not yet supported if the second shape has dimension > 1")
# We cannot make a new functional.PSNEval in
# the natural way, since it tries to compute
# the normal vector by itself.
# Instead, we create things manually, and
# call Functional() with these arguments
sd = ref_el.get_spatial_dimension()
# The pt_dict is a one-item dictionary containing
# the details of the functional.
# The key is the spatial coordinate, which
# is just a concatenation of the two parts.
# The value is a list of tuples, representing
# the normal vector (scaled by the volume of
# the facet) at that point.
# Each tuple looks like (foo, (i,)); the i'th
# component of the scaled normal is foo.
# The following line is only valid when the second
# shape has spatial dimension 1 (enforced above)
Apoint, Avalue = _first_point_pair(Anode)
pt_dict = {Apoint + _first_point(Bnode): Avalue + [(0.0, (len(Apoint),))]}
# The following line should be used in the
# general case
# pt_dict = {Anode.get_point_dict().keys()[0] + Bnode.get_point_dict().keys()[0]: Anode.get_point_dict().values()[0] + [(0.0, (ii,)) for ii in range(len(Anode.get_point_dict().keys()[0]), len(Anode.get_point_dict().keys()[0]) + len(Bnode.get_point_dict().keys()[0]))]}
# THE FOLLOWING IS PROBABLY CORRECT BUT UNTESTED
shp = (sd,)
nodes.append(functional.Functional(ref_el, shp, pt_dict, {}, "PointScaledNormalEval"))
else:
raise NotImplementedError("unsupported functional type")
elif isinstance(Anode, functional.PointEdgeTangentEvaluation):
for Bnode in Bnodes:
if isinstance(Bnode, functional.PointEvaluation):
# case: PointEdgeTangentEval x PointEval
# this is very similar to the case above, so comments omitted
if len(_first_point(Bnode)) > 1:
raise NotImplementedError("PointEdgeTangentEval x PointEval is not yet supported if the second shape has dimension > 1")
sd = ref_el.get_spatial_dimension()
Apoint, Avalue = _first_point_pair(Anode)
pt_dict = {Apoint + _first_point(Bnode): Avalue + [(0.0, (len(Apoint),))]}
# THE FOLLOWING IS PROBABLY CORRECT BUT UNTESTED
shp = (sd,)
nodes.append(functional.Functional(ref_el, shp, pt_dict, {}, "PointEdgeTangent"))
else:
raise NotImplementedError("unsupported functional type")
elif isinstance(Anode, functional.ComponentPointEvaluation):
for Bnode in Bnodes:
if isinstance(Bnode, functional.PointEvaluation):
# case: ComponentPointEval x PointEval
# the CptPointEval functional requires the component
# and the coordinates. very similar to PE x PE case.
sd = ref_el.get_spatial_dimension()
nodes.append(functional.ComponentPointEvaluation(ref_el, Anode.comp, (sd,), _first_point(Anode) + _first_point(Bnode)))
else:
raise NotImplementedError("unsupported functional type")
elif isinstance(Anode, functional.FrobeniusIntegralMoment):
for Bnode in Bnodes:
if isinstance(Bnode, functional.PointEvaluation):
# case: FroIntMom x PointEval
sd = ref_el.get_spatial_dimension()
pt_dict = {}
pt_old = Anode.get_point_dict()
for pt in pt_old:
pt_dict[pt+_first_point(Bnode)] = pt_old[pt] + [(0.0, sd-1)]
# THE FOLLOWING IS PROBABLY CORRECT BUT UNTESTED
shp = (sd,)
nodes.append(functional.Functional(ref_el, shp, pt_dict, {}, "FrobeniusIntegralMoment"))
else:
raise NotImplementedError("unsupported functional type")
elif isinstance(Anode, functional.IntegralMoment):
for Bnode in Bnodes:
if isinstance(Bnode, functional.PointEvaluation):
# case: IntMom x PointEval
sd = ref_el.get_spatial_dimension()
pt_dict = {}
pt_old = Anode.get_point_dict()
for pt in pt_old:
pt_dict[pt+_first_point(Bnode)] = pt_old[pt]
# THE FOLLOWING IS PROBABLY CORRECT BUT UNTESTED
shp = (sd,)
nodes.append(functional.Functional(ref_el, shp, pt_dict, {}, "IntegralMoment"))
else:
> raise NotImplementedError("unsupported functional type")
E NotImplementedError: unsupported functional type
../firedrake_venv/src/fiat/FIAT/tensor_product.py:198: NotImplementedError
Check failure on line 329 in tests/regression/test_fdm.py
github-actions / Firedrake real
test_fdm.test_ipdg_direct_solver[cg-Rectangle]
NotImplementedError: unsupported functional type
Raw output
fs = WithGeometry(FunctionSpace(<firedrake.mesh.MeshTopology object at 0x7fdfa82e7da0>, VectorElement(FiniteElement('Q', qu...ipdg'), dim=3, variant='fdm_ipdg'), name=None), Mesh(VectorElement(FiniteElement('Q', quadrilateral, 1), dim=2), 3982))
@pytest.mark.skipcomplex
def test_ipdg_direct_solver(fs):
mesh = fs.mesh()
x = SpatialCoordinate(mesh)
gdim = mesh.geometric_dimension()
ncomp = fs.ufl_element().value_size
homogenize = gdim > 2
if homogenize:
rg = RandomGenerator(PCG64(seed=123456789))
uh = rg.uniform(fs, -1, 1)
u_exact = zero(uh.ufl_shape)
u_bc = 0
else:
uh = Function(fs)
u_exact = dot(x, x)
if ncomp:
u_exact = as_vector([u_exact + Constant(k) for k in range(ncomp)])
u_bc = u_exact
degree = max(as_tuple(fs.ufl_element().degree()))
quad_degree = 2*(degree+1)-1
# problem coefficients
A1 = diag(Constant(range(1, gdim+1)))
A2 = diag(Constant(range(1, ncomp+1)))
alpha = lambda grad_u: dot(dot(A2, grad_u), A1)
beta = diag(Constant(range(2, ncomp+2)))
extruded = mesh.cell_set._extruded
subs = (1,)
if gdim > 1:
subs += (3,)
if extruded:
subs += ("top",)
bcs = [DirichletBC(fs, u_bc, sub) for sub in subs]
dirichlet_ids = subs
if "on_boundary" in dirichlet_ids:
neumann_ids = []
else:
make_tuple = lambda s: s if type(s) == tuple else (s,)
neumann_ids = list(set(mesh.exterior_facets.unique_markers) - set(sum([make_tuple(s) for s in subs if type(s) != str], ())))
if extruded:
if "top" not in dirichlet_ids:
neumann_ids.append("top")
if "bottom" not in dirichlet_ids:
neumann_ids.append("bottom")
dxq = dx(degree=quad_degree, domain=mesh)
if extruded:
dS_int = dS_v(degree=quad_degree) + dS_h(degree=quad_degree)
ds_ext = {"on_boundary": ds_v(degree=quad_degree), "bottom": ds_b(degree=quad_degree), "top": ds_t(degree=quad_degree)}
ds_Dir = [ds_ext.get(s) or ds_v(s, degree=quad_degree) for s in dirichlet_ids]
ds_Neu = [ds_ext.get(s) or ds_v(s, degree=quad_degree) for s in neumann_ids]
else:
dS_int = dS(degree=quad_degree)
ds_ext = {"on_boundary": ds(degree=quad_degree)}
ds_Dir = [ds_ext.get(s) or ds(s, degree=quad_degree) for s in dirichlet_ids]
ds_Neu = [ds_ext.get(s) or ds(s, degree=quad_degree) for s in neumann_ids]
ds_Dir = sum(ds_Dir, ds(tuple()))
ds_Neu = sum(ds_Neu, ds(tuple()))
n = FacetNormal(mesh)
h = CellVolume(mesh) / FacetArea(mesh)
eta = Constant((degree+1)**2)
penalty = eta / h
num_flux = lambda u: avg(penalty) * avg(outer(u, n))
num_flux_b = lambda u: (penalty/2) * outer(u, n)
a_int = lambda v, u: inner(2 * avg(outer(v, n)), alpha(num_flux(u) - avg(grad(u))))
a_Dir = lambda v, u: inner(outer(v, n), alpha(num_flux_b(u) - grad(u)))
v = TestFunction(fs)
u = TrialFunction(fs)
a = ((inner(v, dot(beta, u)) + inner(grad(v), alpha(grad(u)))) * dxq
+ (a_int(v, u) + a_int(u, v)) * dS_int
+ (a_Dir(v, u) + a_Dir(u, v)) * ds_Dir)
if homogenize:
L = 0
else:
f_exact = alpha(grad(u_exact))
B = dot(beta, u_exact) - div(f_exact)
T = dot(f_exact, n)
L = (inner(v, B)*dxq + inner(v, T)*ds_Neu
+ inner(outer(u_exact, n), alpha(2*num_flux_b(v) - grad(v))) * ds_Dir)
problem = LinearVariationalProblem(a, L, uh, bcs=bcs)
solver = LinearVariationalSolver(problem, solver_parameters={
"mat_type": "matfree",
"ksp_type": "cg",
"ksp_atol": 0.0E0,
"ksp_rtol": 1.0E-8,
"ksp_max_it": 3,
"ksp_monitor": None,
"ksp_norm_type": "unpreconditioned",
"pc_type": "python",
"pc_python_type": "firedrake.PoissonFDMPC",
"fdm_pc_type": "cholesky",
"fdm_pc_factor_mat_solver_type": DEFAULT_DIRECT_SOLVER,
"fdm_pc_factor_mat_ordering_type": "nd",
}, appctx={"eta": eta, })
> solver.solve()
tests/regression/test_fdm.py:329:
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
petsc4py/PETSc/Log.pyx:188: in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
???
petsc4py/PETSc/Log.pyx:189: in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
???
firedrake/adjoint_utils/variational_solver.py:89: in wrapper
out = solve(self, **kwargs)
firedrake/variational_solver.py:307: in solve
dbc.apply(problem.u_restrict)
petsc4py/PETSc/Log.pyx:188: in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
???
petsc4py/PETSc/Log.pyx:189: in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
???
firedrake/adjoint_utils/dirichletbc.py:32: in wrapper
ret = apply(self, *args, **kwargs)
firedrake/bcs.py:451: in apply
r.assign(self.function_arg, subset=self.node_set)
../firedrake_venv/src/PyOP2/pyop2/utils.py:62: in __get__
obj.__dict__[self.__name__] = result = self.fget(obj)
firedrake/bcs.py:179: in node_set
return op2.Subset(self._function_space.node_set, self.nodes)
../firedrake_venv/src/PyOP2/pyop2/utils.py:62: in __get__
obj.__dict__[self.__name__] = result = self.fget(obj)
firedrake/bcs.py:169: in nodes
bcnodes1.append(hermite_stride(self._function_space.boundary_nodes(ss)))
firedrake/bcs.py:143: in hermite_stride
deriv_nodes = [i for i, node in enumerate(fe.fiat_equivalent.dual) if len(node.deriv_dict) != 0]
../firedrake_venv/src/tsfc/gem/utils.py:19: in __get__
obj.__dict__[self.__name__] = result = self.fget(obj)
../firedrake_venv/src/FInAT/finat/cube.py:65: in fiat_equivalent
return FIAT_FlattenedDimensions(self.product.fiat_equivalent)
../firedrake_venv/src/tsfc/gem/utils.py:19: in __get__
obj.__dict__[self.__name__] = result = self.fget(obj)
../firedrake_venv/src/FInAT/finat/tensor_product.py:82: in fiat_equivalent
return FIAT.TensorProductElement(A.fiat_equivalent, B.fiat_equivalent)
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
self = <FIAT.tensor_product.TensorProductElement object at 0x7fdfa839c710>
A = <FIAT.fdm_element.FDMLagrange object at 0x7fdfa82e4410>
B = <FIAT.fdm_element.FDMLagrange object at 0x7fdfa82e4410>
def __init__(self, A, B):
# set up simple things
order = min(A.get_order(), B.get_order())
if A.get_formdegree() is None or B.get_formdegree() is None:
formdegree = None
else:
formdegree = A.get_formdegree() + B.get_formdegree()
# set up reference element
ref_el = TensorProductCell(A.get_reference_element(),
B.get_reference_element())
if A.mapping()[0] != "affine" and B.mapping()[0] == "affine":
mapping = A.mapping()[0]
elif B.mapping()[0] != "affine" and A.mapping()[0] == "affine":
mapping = B.mapping()[0]
elif A.mapping()[0] == "affine" and B.mapping()[0] == "affine":
mapping = "affine"
else:
raise ValueError("check tensor product mappings - at least one must be affine")
# set up entity_ids
Adofs = A.entity_dofs()
Bdofs = B.entity_dofs()
Bsdim = B.space_dimension()
entity_ids = {}
for curAdim in Adofs:
for curBdim in Bdofs:
entity_ids[(curAdim, curBdim)] = {}
dim_cur = 0
for entityA in Adofs[curAdim]:
for entityB in Bdofs[curBdim]:
entity_ids[(curAdim, curBdim)][dim_cur] = \
[x*Bsdim + y for x in Adofs[curAdim][entityA]
for y in Bdofs[curBdim][entityB]]
dim_cur += 1
# set up dual basis
Anodes = A.dual_basis()
Bnodes = B.dual_basis()
# build the dual set by inspecting the current dual
# sets item by item.
# Currently supported cases:
# PointEval x PointEval = PointEval [scalar x scalar = scalar]
# PointScaledNormalEval x PointEval = PointScaledNormalEval [vector x scalar = vector]
# ComponentPointEvaluation x PointEval [vector x scalar = vector]
nodes = []
for Anode in Anodes:
if isinstance(Anode, functional.PointEvaluation):
for Bnode in Bnodes:
if isinstance(Bnode, functional.PointEvaluation):
# case: PointEval x PointEval
# the PointEval functional just requires the
# coordinates. these are currently stored as
# the key of a one-item dictionary. we retrieve
# these by calling get_point_dict(), and
# use the concatenation to make a new PointEval
nodes.append(functional.PointEvaluation(ref_el, _first_point(Anode) + _first_point(Bnode)))
elif isinstance(Bnode, functional.IntegralMoment):
# dummy functional for product with integral moments
nodes.append(functional.Functional(None, None, None,
{}, "Undefined"))
elif isinstance(Bnode, functional.PointDerivative):
# dummy functional for product with point derivative
nodes.append(functional.Functional(None, None, None,
{}, "Undefined"))
else:
raise NotImplementedError("unsupported functional type")
elif isinstance(Anode, functional.PointScaledNormalEvaluation):
for Bnode in Bnodes:
if isinstance(Bnode, functional.PointEvaluation):
# case: PointScaledNormalEval x PointEval
# this could be wrong if the second shape
# has spatial dimension >1, since we are not
# explicitly scaling by facet size
if len(_first_point(Bnode)) > 1:
# TODO: support this case one day
raise NotImplementedError("PointScaledNormalEval x PointEval is not yet supported if the second shape has dimension > 1")
# We cannot make a new functional.PSNEval in
# the natural way, since it tries to compute
# the normal vector by itself.
# Instead, we create things manually, and
# call Functional() with these arguments
sd = ref_el.get_spatial_dimension()
# The pt_dict is a one-item dictionary containing
# the details of the functional.
# The key is the spatial coordinate, which
# is just a concatenation of the two parts.
# The value is a list of tuples, representing
# the normal vector (scaled by the volume of
# the facet) at that point.
# Each tuple looks like (foo, (i,)); the i'th
# component of the scaled normal is foo.
# The following line is only valid when the second
# shape has spatial dimension 1 (enforced above)
Apoint, Avalue = _first_point_pair(Anode)
pt_dict = {Apoint + _first_point(Bnode): Avalue + [(0.0, (len(Apoint),))]}
# The following line should be used in the
# general case
# pt_dict = {Anode.get_point_dict().keys()[0] + Bnode.get_point_dict().keys()[0]: Anode.get_point_dict().values()[0] + [(0.0, (ii,)) for ii in range(len(Anode.get_point_dict().keys()[0]), len(Anode.get_point_dict().keys()[0]) + len(Bnode.get_point_dict().keys()[0]))]}
# THE FOLLOWING IS PROBABLY CORRECT BUT UNTESTED
shp = (sd,)
nodes.append(functional.Functional(ref_el, shp, pt_dict, {}, "PointScaledNormalEval"))
else:
raise NotImplementedError("unsupported functional type")
elif isinstance(Anode, functional.PointEdgeTangentEvaluation):
for Bnode in Bnodes:
if isinstance(Bnode, functional.PointEvaluation):
# case: PointEdgeTangentEval x PointEval
# this is very similar to the case above, so comments omitted
if len(_first_point(Bnode)) > 1:
raise NotImplementedError("PointEdgeTangentEval x PointEval is not yet supported if the second shape has dimension > 1")
sd = ref_el.get_spatial_dimension()
Apoint, Avalue = _first_point_pair(Anode)
pt_dict = {Apoint + _first_point(Bnode): Avalue + [(0.0, (len(Apoint),))]}
# THE FOLLOWING IS PROBABLY CORRECT BUT UNTESTED
shp = (sd,)
nodes.append(functional.Functional(ref_el, shp, pt_dict, {}, "PointEdgeTangent"))
else:
raise NotImplementedError("unsupported functional type")
elif isinstance(Anode, functional.ComponentPointEvaluation):
for Bnode in Bnodes:
if isinstance(Bnode, functional.PointEvaluation):
# case: ComponentPointEval x PointEval
# the CptPointEval functional requires the component
# and the coordinates. very similar to PE x PE case.
sd = ref_el.get_spatial_dimension()
nodes.append(functional.ComponentPointEvaluation(ref_el, Anode.comp, (sd,), _first_point(Anode) + _first_point(Bnode)))
else:
raise NotImplementedError("unsupported functional type")
elif isinstance(Anode, functional.FrobeniusIntegralMoment):
for Bnode in Bnodes:
if isinstance(Bnode, functional.PointEvaluation):
# case: FroIntMom x PointEval
sd = ref_el.get_spatial_dimension()
pt_dict = {}
pt_old = Anode.get_point_dict()
for pt in pt_old:
pt_dict[pt+_first_point(Bnode)] = pt_old[pt] + [(0.0, sd-1)]
# THE FOLLOWING IS PROBABLY CORRECT BUT UNTESTED
shp = (sd,)
nodes.append(functional.Functional(ref_el, shp, pt_dict, {}, "FrobeniusIntegralMoment"))
else:
raise NotImplementedError("unsupported functional type")
elif isinstance(Anode, functional.IntegralMoment):
for Bnode in Bnodes:
if isinstance(Bnode, functional.PointEvaluation):
# case: IntMom x PointEval
sd = ref_el.get_spatial_dimension()
pt_dict = {}
pt_old = Anode.get_point_dict()
for pt in pt_old:
pt_dict[pt+_first_point(Bnode)] = pt_old[pt]
# THE FOLLOWING IS PROBABLY CORRECT BUT UNTESTED
shp = (sd,)
nodes.append(functional.Functional(ref_el, shp, pt_dict, {}, "IntegralMoment"))
else:
> raise NotImplementedError("unsupported functional type")
E NotImplementedError: unsupported functional type
../firedrake_venv/src/fiat/FIAT/tensor_product.py:198: NotImplementedError
Check failure on line 329 in tests/regression/test_fdm.py
github-actions / Firedrake real
test_fdm.test_ipdg_direct_solver[cg-Box]
NotImplementedError: unsupported functional type
Raw output
fs = WithGeometry(FunctionSpace(<firedrake.mesh.ExtrudedMeshTopology object at 0x7fdfa8b7ec60>, VectorElement(TensorProduct...rilateral, 1), FiniteElement('Lagrange', interval, 1), cell=TensorProductCell(quadrilateral, interval)), dim=3), 3999))
@pytest.mark.skipcomplex
def test_ipdg_direct_solver(fs):
mesh = fs.mesh()
x = SpatialCoordinate(mesh)
gdim = mesh.geometric_dimension()
ncomp = fs.ufl_element().value_size
homogenize = gdim > 2
if homogenize:
rg = RandomGenerator(PCG64(seed=123456789))
uh = rg.uniform(fs, -1, 1)
u_exact = zero(uh.ufl_shape)
u_bc = 0
else:
uh = Function(fs)
u_exact = dot(x, x)
if ncomp:
u_exact = as_vector([u_exact + Constant(k) for k in range(ncomp)])
u_bc = u_exact
degree = max(as_tuple(fs.ufl_element().degree()))
quad_degree = 2*(degree+1)-1
# problem coefficients
A1 = diag(Constant(range(1, gdim+1)))
A2 = diag(Constant(range(1, ncomp+1)))
alpha = lambda grad_u: dot(dot(A2, grad_u), A1)
beta = diag(Constant(range(2, ncomp+2)))
extruded = mesh.cell_set._extruded
subs = (1,)
if gdim > 1:
subs += (3,)
if extruded:
subs += ("top",)
bcs = [DirichletBC(fs, u_bc, sub) for sub in subs]
dirichlet_ids = subs
if "on_boundary" in dirichlet_ids:
neumann_ids = []
else:
make_tuple = lambda s: s if type(s) == tuple else (s,)
neumann_ids = list(set(mesh.exterior_facets.unique_markers) - set(sum([make_tuple(s) for s in subs if type(s) != str], ())))
if extruded:
if "top" not in dirichlet_ids:
neumann_ids.append("top")
if "bottom" not in dirichlet_ids:
neumann_ids.append("bottom")
dxq = dx(degree=quad_degree, domain=mesh)
if extruded:
dS_int = dS_v(degree=quad_degree) + dS_h(degree=quad_degree)
ds_ext = {"on_boundary": ds_v(degree=quad_degree), "bottom": ds_b(degree=quad_degree), "top": ds_t(degree=quad_degree)}
ds_Dir = [ds_ext.get(s) or ds_v(s, degree=quad_degree) for s in dirichlet_ids]
ds_Neu = [ds_ext.get(s) or ds_v(s, degree=quad_degree) for s in neumann_ids]
else:
dS_int = dS(degree=quad_degree)
ds_ext = {"on_boundary": ds(degree=quad_degree)}
ds_Dir = [ds_ext.get(s) or ds(s, degree=quad_degree) for s in dirichlet_ids]
ds_Neu = [ds_ext.get(s) or ds(s, degree=quad_degree) for s in neumann_ids]
ds_Dir = sum(ds_Dir, ds(tuple()))
ds_Neu = sum(ds_Neu, ds(tuple()))
n = FacetNormal(mesh)
h = CellVolume(mesh) / FacetArea(mesh)
eta = Constant((degree+1)**2)
penalty = eta / h
num_flux = lambda u: avg(penalty) * avg(outer(u, n))
num_flux_b = lambda u: (penalty/2) * outer(u, n)
a_int = lambda v, u: inner(2 * avg(outer(v, n)), alpha(num_flux(u) - avg(grad(u))))
a_Dir = lambda v, u: inner(outer(v, n), alpha(num_flux_b(u) - grad(u)))
v = TestFunction(fs)
u = TrialFunction(fs)
a = ((inner(v, dot(beta, u)) + inner(grad(v), alpha(grad(u)))) * dxq
+ (a_int(v, u) + a_int(u, v)) * dS_int
+ (a_Dir(v, u) + a_Dir(u, v)) * ds_Dir)
if homogenize:
L = 0
else:
f_exact = alpha(grad(u_exact))
B = dot(beta, u_exact) - div(f_exact)
T = dot(f_exact, n)
L = (inner(v, B)*dxq + inner(v, T)*ds_Neu
+ inner(outer(u_exact, n), alpha(2*num_flux_b(v) - grad(v))) * ds_Dir)
problem = LinearVariationalProblem(a, L, uh, bcs=bcs)
solver = LinearVariationalSolver(problem, solver_parameters={
"mat_type": "matfree",
"ksp_type": "cg",
"ksp_atol": 0.0E0,
"ksp_rtol": 1.0E-8,
"ksp_max_it": 3,
"ksp_monitor": None,
"ksp_norm_type": "unpreconditioned",
"pc_type": "python",
"pc_python_type": "firedrake.PoissonFDMPC",
"fdm_pc_type": "cholesky",
"fdm_pc_factor_mat_solver_type": DEFAULT_DIRECT_SOLVER,
"fdm_pc_factor_mat_ordering_type": "nd",
}, appctx={"eta": eta, })
> solver.solve()
tests/regression/test_fdm.py:329:
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
petsc4py/PETSc/Log.pyx:188: in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
???
petsc4py/PETSc/Log.pyx:189: in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
???
firedrake/adjoint_utils/variational_solver.py:89: in wrapper
out = solve(self, **kwargs)
firedrake/variational_solver.py:307: in solve
dbc.apply(problem.u_restrict)
petsc4py/PETSc/Log.pyx:188: in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
???
petsc4py/PETSc/Log.pyx:189: in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
???
firedrake/adjoint_utils/dirichletbc.py:32: in wrapper
ret = apply(self, *args, **kwargs)
firedrake/bcs.py:451: in apply
r.assign(self.function_arg, subset=self.node_set)
../firedrake_venv/src/PyOP2/pyop2/utils.py:62: in __get__
obj.__dict__[self.__name__] = result = self.fget(obj)
firedrake/bcs.py:179: in node_set
return op2.Subset(self._function_space.node_set, self.nodes)
../firedrake_venv/src/PyOP2/pyop2/utils.py:62: in __get__
obj.__dict__[self.__name__] = result = self.fget(obj)
firedrake/bcs.py:169: in nodes
bcnodes1.append(hermite_stride(self._function_space.boundary_nodes(ss)))
firedrake/bcs.py:143: in hermite_stride
deriv_nodes = [i for i, node in enumerate(fe.fiat_equivalent.dual) if len(node.deriv_dict) != 0]
../firedrake_venv/src/tsfc/gem/utils.py:19: in __get__
obj.__dict__[self.__name__] = result = self.fget(obj)
../firedrake_venv/src/FInAT/finat/tensor_product.py:82: in fiat_equivalent
return FIAT.TensorProductElement(A.fiat_equivalent, B.fiat_equivalent)
../firedrake_venv/src/tsfc/gem/utils.py:19: in __get__
obj.__dict__[self.__name__] = result = self.fget(obj)
../firedrake_venv/src/FInAT/finat/cube.py:65: in fiat_equivalent
return FIAT_FlattenedDimensions(self.product.fiat_equivalent)
../firedrake_venv/src/tsfc/gem/utils.py:19: in __get__
obj.__dict__[self.__name__] = result = self.fget(obj)
../firedrake_venv/src/FInAT/finat/tensor_product.py:82: in fiat_equivalent
return FIAT.TensorProductElement(A.fiat_equivalent, B.fiat_equivalent)
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
self = <FIAT.tensor_product.TensorProductElement object at 0x7fdfa7425340>
A = <FIAT.fdm_element.FDMLagrange object at 0x7fdfa82e4410>
B = <FIAT.fdm_element.FDMLagrange object at 0x7fdfa82e4410>
def __init__(self, A, B):
# set up simple things
order = min(A.get_order(), B.get_order())
if A.get_formdegree() is None or B.get_formdegree() is None:
formdegree = None
else:
formdegree = A.get_formdegree() + B.get_formdegree()
# set up reference element
ref_el = TensorProductCell(A.get_reference_element(),
B.get_reference_element())
if A.mapping()[0] != "affine" and B.mapping()[0] == "affine":
mapping = A.mapping()[0]
elif B.mapping()[0] != "affine" and A.mapping()[0] == "affine":
mapping = B.mapping()[0]
elif A.mapping()[0] == "affine" and B.mapping()[0] == "affine":
mapping = "affine"
else:
raise ValueError("check tensor product mappings - at least one must be affine")
# set up entity_ids
Adofs = A.entity_dofs()
Bdofs = B.entity_dofs()
Bsdim = B.space_dimension()
entity_ids = {}
for curAdim in Adofs:
for curBdim in Bdofs:
entity_ids[(curAdim, curBdim)] = {}
dim_cur = 0
for entityA in Adofs[curAdim]:
for entityB in Bdofs[curBdim]:
entity_ids[(curAdim, curBdim)][dim_cur] = \
[x*Bsdim + y for x in Adofs[curAdim][entityA]
for y in Bdofs[curBdim][entityB]]
dim_cur += 1
# set up dual basis
Anodes = A.dual_basis()
Bnodes = B.dual_basis()
# build the dual set by inspecting the current dual
# sets item by item.
# Currently supported cases:
# PointEval x PointEval = PointEval [scalar x scalar = scalar]
# PointScaledNormalEval x PointEval = PointScaledNormalEval [vector x scalar = vector]
# ComponentPointEvaluation x PointEval [vector x scalar = vector]
nodes = []
for Anode in Anodes:
if isinstance(Anode, functional.PointEvaluation):
for Bnode in Bnodes:
if isinstance(Bnode, functional.PointEvaluation):
# case: PointEval x PointEval
# the PointEval functional just requires the
# coordinates. these are currently stored as
# the key of a one-item dictionary. we retrieve
# these by calling get_point_dict(), and
# use the concatenation to make a new PointEval
nodes.append(functional.PointEvaluation(ref_el, _first_point(Anode) + _first_point(Bnode)))
elif isinstance(Bnode, functional.IntegralMoment):
# dummy functional for product with integral moments
nodes.append(functional.Functional(None, None, None,
{}, "Undefined"))
elif isinstance(Bnode, functional.PointDerivative):
# dummy functional for product with point derivative
nodes.append(functional.Functional(None, None, None,
{}, "Undefined"))
else:
raise NotImplementedError("unsupported functional type")
elif isinstance(Anode, functional.PointScaledNormalEvaluation):
for Bnode in Bnodes:
if isinstance(Bnode, functional.PointEvaluation):
# case: PointScaledNormalEval x PointEval
# this could be wrong if the second shape
# has spatial dimension >1, since we are not
# explicitly scaling by facet size
if len(_first_point(Bnode)) > 1:
# TODO: support this case one day
raise NotImplementedError("PointScaledNormalEval x PointEval is not yet supported if the second shape has dimension > 1")
# We cannot make a new functional.PSNEval in
# the natural way, since it tries to compute
# the normal vector by itself.
# Instead, we create things manually, and
# call Functional() with these arguments
sd = ref_el.get_spatial_dimension()
# The pt_dict is a one-item dictionary containing
# the details of the functional.
# The key is the spatial coordinate, which
# is just a concatenation of the two parts.
# The value is a list of tuples, representing
# the normal vector (scaled by the volume of
# the facet) at that point.
# Each tuple looks like (foo, (i,)); the i'th
# component of the scaled normal is foo.
# The following line is only valid when the second
# shape has spatial dimension 1 (enforced above)
Apoint, Avalue = _first_point_pair(Anode)
pt_dict = {Apoint + _first_point(Bnode): Avalue + [(0.0, (len(Apoint),))]}
# The following line should be used in the
# general case
# pt_dict = {Anode.get_point_dict().keys()[0] + Bnode.get_point_dict().keys()[0]: Anode.get_point_dict().values()[0] + [(0.0, (ii,)) for ii in range(len(Anode.get_point_dict().keys()[0]), len(Anode.get_point_dict().keys()[0]) + len(Bnode.get_point_dict().keys()[0]))]}
# THE FOLLOWING IS PROBABLY CORRECT BUT UNTESTED
shp = (sd,)
nodes.append(functional.Functional(ref_el, shp, pt_dict, {}, "PointScaledNormalEval"))
else:
raise NotImplementedError("unsupported functional type")
elif isinstance(Anode, functional.PointEdgeTangentEvaluation):
for Bnode in Bnodes:
if isinstance(Bnode, functional.PointEvaluation):
# case: PointEdgeTangentEval x PointEval
# this is very similar to the case above, so comments omitted
if len(_first_point(Bnode)) > 1:
raise NotImplementedError("PointEdgeTangentEval x PointEval is not yet supported if the second shape has dimension > 1")
sd = ref_el.get_spatial_dimension()
Apoint, Avalue = _first_point_pair(Anode)
pt_dict = {Apoint + _first_point(Bnode): Avalue + [(0.0, (len(Apoint),))]}
# THE FOLLOWING IS PROBABLY CORRECT BUT UNTESTED
shp = (sd,)
nodes.append(functional.Functional(ref_el, shp, pt_dict, {}, "PointEdgeTangent"))
else:
raise NotImplementedError("unsupported functional type")
elif isinstance(Anode, functional.ComponentPointEvaluation):
for Bnode in Bnodes:
if isinstance(Bnode, functional.PointEvaluation):
# case: ComponentPointEval x PointEval
# the CptPointEval functional requires the component
# and the coordinates. very similar to PE x PE case.
sd = ref_el.get_spatial_dimension()
nodes.append(functional.ComponentPointEvaluation(ref_el, Anode.comp, (sd,), _first_point(Anode) + _first_point(Bnode)))
else:
raise NotImplementedError("unsupported functional type")
elif isinstance(Anode, functional.FrobeniusIntegralMoment):
for Bnode in Bnodes:
if isinstance(Bnode, functional.PointEvaluation):
# case: FroIntMom x PointEval
sd = ref_el.get_spatial_dimension()
pt_dict = {}
pt_old = Anode.get_point_dict()
for pt in pt_old:
pt_dict[pt+_first_point(Bnode)] = pt_old[pt] + [(0.0, sd-1)]
# THE FOLLOWING IS PROBABLY CORRECT BUT UNTESTED
shp = (sd,)
nodes.append(functional.Functional(ref_el, shp, pt_dict, {}, "FrobeniusIntegralMoment"))
else:
raise NotImplementedError("unsupported functional type")
elif isinstance(Anode, functional.IntegralMoment):
for Bnode in Bnodes:
if isinstance(Bnode, functional.PointEvaluation):
# case: IntMom x PointEval
sd = ref_el.get_spatial_dimension()
pt_dict = {}
pt_old = Anode.get_point_dict()
for pt in pt_old:
pt_dict[pt+_first_point(Bnode)] = pt_old[pt]
# THE FOLLOWING IS PROBABLY CORRECT BUT UNTESTED
shp = (sd,)
nodes.append(functional.Functional(ref_el, shp, pt_dict, {}, "IntegralMoment"))
else:
> raise NotImplementedError("unsupported functional type")
E NotImplementedError: unsupported functional type
../firedrake_venv/src/fiat/FIAT/tensor_product.py:198: NotImplementedError