The discretization of the following BilinearForm, when domain is a mapped domain, does not work.
q_half = BilinearForm((u1, v1), integral(domain, dot( u1, cross(F, v1))))
Notably, the discretization of the similar looking BilinearForm
q = BilinearForm((u1, v1), integral(domain, dot( cross(F, u1), cross(F, v1) )))
does work. This issue is not Pyccel related, as it occurs both with and without passing a backend variable.
I need to be able to discretize q_half to (hopefully) improve the convergence behaviour of my MHD code.
The Traceback error message is rather long, and from my understanding of it, the issue is sympy related.
I include both a MWE and the corresponding error message.
import time
from mpi4py import MPI
import numpy as np
from sympde.calculus import dot, cross
from sympde.expr import BilinearForm, integral
from sympde.topology import elements_of, Cube, Mapping, Derham
from psydac.api.discretization import discretize
from psydac.api.settings import PSYDAC_BACKEND_GPYCCEL
from psydac.fem.basic import FemField
comm = MPI.COMM_WORLD
backend = PSYDAC_BACKEND_GPYCCEL
# ----- Parameters to play around with -----
ncells = [12, 8, 16]
degree = [2, 4, 3]
periodic = [False, True, False]
mult = [1, 3, 2]
# Choose
# None for no mapping, or (no issue here)
# 'analytical' for an analytical mapping ( issue here)
mapping_option = 'analytical' # None
# ------------------------------------------
#
#
#
logical_domain = Cube('C', bounds1=(0,1), bounds2=(0,1), bounds3=(0,1))
if mapping_option is None:
domain = logical_domain
elif mapping_option == 'analytical':
class SquareTorus(Mapping):
_expressions = {'x': 'x1 * cos(x2)',
'y': 'x1 * sin(x2)',
'z': 'x3'}
_ldim = 3
_pdim = 3
mapping = SquareTorus('ST')
domain = mapping(logical_domain)
else:
raise ValueError(f'mapping_option {mapping_option} not understood.')
derham = Derham(domain)
domain_h = discretize(domain, ncells=ncells, periodic=periodic, comm=comm)
derham_h = discretize(derham, domain_h, degree=degree, multiplicity=mult)
V1 = derham.V1
V1h = derham_h.V1
u1, v1, F = elements_of(V1, names='u1, v1, F')
# generate a random free FemField
F_coeffs = V1h.coeff_space.zeros()
rng = np.random.default_rng(seed=42)
for block in F_coeffs.blocks:
rng.random(size=block._data.shape, dtype="float64", out=block._data)
F_field = FemField(V1h, F_coeffs)
# Previously I worked with q in my code. A version with (hopefully) better convergence behaviour requires q_half
q = BilinearForm((u1, v1), integral(domain, dot( cross(F, u1), cross(F, v1) )))
q_half = BilinearForm((u1, v1), integral(domain, dot( u1, cross(F, v1))))
# Discretizing q works just fine
t0 = time.time()
q_h = discretize(q, domain_h, (V1h, V1h), backend=backend)
t1 = time.time()
print(f'q discretized in {t1-t0:.3g}')
print(f' ----- ')
print()
# Discretizing q_half does not work, unless not using a mapping, and the backend is not the problem!
t0 = time.time()
q_half_h = discretize(q_half, domain_h, (V1h, V1h), backend=backend)
t1 = time.time()
print(f'q_half discretized in {t1-t0:.3g}')
print(f' ----- ')
print()Executing this code results in the following error message:
q discretized in 18.6
-----
Traceback (most recent call last):
File "/home/jubuntu/psydac-hodge/.venv/lib/python3.11/site-packages/sympy/polys/matrices/sdm.py", line 80, in getitem
return self[i][j]
~~~~~~~^^^
KeyError: 1
During handling of the above exception, another exception occurred:
Traceback (most recent call last):
File "/home/jubuntu/psydac-hodge/.venv/lib/python3.11/site-packages/sympy/matrices/repmatrix.py", line 700, in _getitem_RepMatrix
return self._rep.getitem_sympy(index_(i), index_(j))
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/home/jubuntu/psydac-hodge/.venv/lib/python3.11/site-packages/sympy/polys/matrices/domainmatrix.py", line 146, in getitem_sympy
return self.domain.to_sympy(self.rep.getitem(i, j))
^^^^^^^^^^^^^^^^^^^^^^
File "/home/jubuntu/psydac-hodge/.venv/lib/python3.11/site-packages/sympy/polys/matrices/sdm.py", line 89, in getitem
raise IndexError("index out of range")
IndexError: index out of range
During handling of the above exception, another exception occurred:
Traceback (most recent call last):
File "/home/jubuntu/psydac/psydac/api/tests/mwe.py", line 87, in <module>
q_half_h = discretize(q_half, domain_h, (V1h, V1h), backend=backend)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/home/jubuntu/psydac/psydac/api/discretization.py", line 607, in discretize
kernel_expr = TerminalExpr(a, domain)
^^^^^^^^^^^^^^^^^^^^^^^
File "/home/jubuntu/psydac-hodge/.venv/lib/python3.11/site-packages/sympde/expr/evaluation.py", line 498, in __new__
r = cls.eval(expr, domain)
^^^^^^^^^^^^^^^^^^^^^^
File "/home/jubuntu/psydac-hodge/.venv/lib/python3.11/site-packages/sympy/core/cache.py", line 72, in wrapper
retval = cfunc(*args, **kwargs)
^^^^^^^^^^^^^^^^^^^^^^
File "/home/jubuntu/psydac-hodge/.venv/lib/python3.11/site-packages/sympde/expr/evaluation.py", line 702, in eval
newexpr = cls.eval(a, domain=domain)
^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/home/jubuntu/psydac-hodge/.venv/lib/python3.11/site-packages/sympy/core/cache.py", line 72, in wrapper
retval = cfunc(*args, **kwargs)
^^^^^^^^^^^^^^^^^^^^^^
File "/home/jubuntu/psydac-hodge/.venv/lib/python3.11/site-packages/sympde/expr/evaluation.py", line 758, in eval
expr = cls.eval(expr._args[0], domain=domain)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/home/jubuntu/psydac-hodge/.venv/lib/python3.11/site-packages/sympy/core/cache.py", line 72, in wrapper
retval = cfunc(*args, **kwargs)
^^^^^^^^^^^^^^^^^^^^^^
File "/home/jubuntu/psydac-hodge/.venv/lib/python3.11/site-packages/sympde/expr/evaluation.py", line 536, in eval
args = [cls.eval(a, domain=domain) for a in expr.args]
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/home/jubuntu/psydac-hodge/.venv/lib/python3.11/site-packages/sympde/expr/evaluation.py", line 536, in <listcomp>
args = [cls.eval(a, domain=domain) for a in expr.args]
^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/home/jubuntu/psydac-hodge/.venv/lib/python3.11/site-packages/sympy/core/cache.py", line 72, in wrapper
retval = cfunc(*args, **kwargs)
^^^^^^^^^^^^^^^^^^^^^^
File "/home/jubuntu/psydac-hodge/.venv/lib/python3.11/site-packages/sympde/expr/evaluation.py", line 787, in eval
return new(*args)
^^^^^^^^^^
File "/home/jubuntu/psydac-hodge/.venv/lib/python3.11/site-packages/sympde/core/algebra.py", line 103, in __new__
r = cls.eval(*args)
^^^^^^^^^^^^^^^
File "/home/jubuntu/psydac-hodge/.venv/lib/python3.11/site-packages/sympde/core/algebra.py", line 181, in eval
return Tuple(u[0,0]*v[0] + u[0,1]*v[1] + u[0,2]*v[2],
~^^^^^
File "/home/jubuntu/psydac-hodge/.venv/lib/python3.11/site-packages/sympy/matrices/repmatrix.py", line 220, in __getitem__
return _getitem_RepMatrix(self, key)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/home/jubuntu/psydac-hodge/.venv/lib/python3.11/site-packages/sympy/matrices/repmatrix.py", line 721, in _getitem_RepMatrix
return self.extract(i, j)
^^^^^^^^^^^^^^^^^^
File "/home/jubuntu/psydac-hodge/.venv/lib/python3.11/site-packages/sympy/matrices/common.py", line 373, in extract
colsList = [a2idx(k, self.cols) for k in colsList]
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/home/jubuntu/psydac-hodge/.venv/lib/python3.11/site-packages/sympy/matrices/common.py", line 373, in <listcomp>
colsList = [a2idx(k, self.cols) for k in colsList]
^^^^^^^^^^^^^^^^^^^
File "/home/jubuntu/psydac-hodge/.venv/lib/python3.11/site-packages/sympy/matrices/common.py", line 3247, in a2idx
raise IndexError("Index out of range: a[%s]" % (j,))
IndexError: Index out of range: a[1]
The discretization of the following
BilinearForm, whendomainis a mapped domain, does not work.Notably, the discretization of the similar looking
BilinearFormdoes work. This issue is not
Pyccelrelated, as it occurs both with and without passing abackendvariable.I need to be able to discretize
q_halfto (hopefully) improve the convergence behaviour of my MHD code.The Traceback error message is rather long, and from my understanding of it, the issue is
sympyrelated.I include both a MWE and the corresponding error message.
Executing this code results in the following error message: