BilinearForm Discretization Issue 2

[jowezarek]

The discretization of the following BilinearForm, when domain is a mapped domain, does not work.

g       = BilinearForm((u1, v1), integral(domain, inner(curl(cross(u1, F)), curl(cross(v1, F)))))

This issue is not Pyccel related, as it occurs both with and without passing a backend variable.

I need to be able to discretize g to (hopefully) improve the convergence behaviour of my MHD code.

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, inner, curl
from    sympde.expr                 import BilinearForm, integral
from    sympde.topology             import element_of, 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
V2  = derham.V2

V1h = derham_h.V1
V2h = derham_h.V2

u1, v1  = elements_of(V1, names='u1, v1')
F       = element_of (V2, name='F')

# generate a random free FemField
F_coeffs = V2h.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(V2h, F_coeffs)

g       = BilinearForm((u1, v1), integral(domain, inner(curl(cross(u1, F)), curl(cross(v1, F)))))

t0 = time.time()
g_h = discretize(g, domain_h, (V1h, V1h), backend=backend)
t1 = time.time()
print(f'g discretized in {t1-t0:.3g}')
print(f' ----- ')
print()

Traceback (most recent call last):
  File "/home/jubuntu/psydac/psydac/api/tests/mwe2.py", line 82, in <module>
    g_h = discretize(g, domain_h, (V1h, V1h), backend=backend)
          ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/jubuntu/psydac/psydac/api/discretization.py", line 604, in discretize
    a       = LogicalExpr (a, domain)
              ^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/jubuntu/psydac-hodge/.venv/lib/python3.11/site-packages/sympde/topology/mapping.py", line 880, in __new__
    r = cls.eval(expr, domain, **options)
        ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/jubuntu/psydac-hodge/.venv/lib/python3.11/site-packages/sympde/topology/mapping.py", line 1233, in eval
    body    = cls.eval(expr.expr, domain)
              ^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/jubuntu/psydac-hodge/.venv/lib/python3.11/site-packages/sympde/topology/mapping.py", line 1226, in eval
    body   = cls.eval(expr.expr, domain)*det
             ^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/jubuntu/psydac-hodge/.venv/lib/python3.11/site-packages/sympde/topology/mapping.py", line 1090, in eval
    args = [cls.eval(arg, domain) for arg in expr.args]
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/jubuntu/psydac-hodge/.venv/lib/python3.11/site-packages/sympde/topology/mapping.py", line 1090, in <listcomp>
    args = [cls.eval(arg, domain) for arg in expr.args]
            ^^^^^^^^^^^^^^^^^^^^^
  File "/home/jubuntu/psydac-hodge/.venv/lib/python3.11/site-packages/sympde/topology/mapping.py", line 1037, in eval
    raise NotImplementedError('TODO')
NotImplementedError: TODO

[FrederikSchnack]

I took a quick deep dive into the workings of SymPDE and came up with a possible solution.

First of all, the error message suggests that the evaluation of the curl with a cross product inside is not implemented for mapped domains, as the corresponding integral transformation is missing. Thus I added the following case to the eval function of LogicalExpr here:

        elif isinstance(expr, curl):
            arg = expr.args[0]

            if isinstance(mapping, InterfaceMapping):
                if isinstance(arg, MinusInterfaceOperator):
                    arg     = arg.args[0]
                    mapping = mapping.minus
                elif isinstance(arg, PlusInterfaceOperator):
                    arg = arg.args[0]
                    mapping = mapping.plus
                else:
                    raise TypeError(arg)


            if isinstance(arg, (VectorFunction)):
                arg = PullBack(arg, mapping)

            elif isinstance(arg, Cross):
                args = [PullBack(a, mapping) for a in arg.args]
                arg  = Cross(*args)

            else:
                arg = cls.eval(arg, domain)

            if isinstance(arg, PullBack) and isinstance(arg.kind, HcurlSpaceType):
                J   = mapping.jacobian
                arg = arg.test

                if isinstance(expr.args[0], (MinusInterfaceOperator, PlusInterfaceOperator)):
                    arg = type(expr.args[0])(arg)
                if expr.is_scalar:
                    return (1/J.det())*curl(arg)

                return (J/J.det())*curl(arg)

            elif isinstance(arg, Cross):

                arg = Cross(args[0].expr, args[1].expr)
                
                J   = mapping.jacobian

                return J * curl( J.inv() * arg)
            else:
                raise NotImplementedError('TODO')

This leads to another error: TypeError: 'ImmutableDenseMatrix' object cannot be interpreted as an integer. Turns out that when evaluating the matrix-vector product between the Jacobian and the resulting cross-product, the multiplication doesn't work. To fix this, I changed the eval function of Cross_3d here to actually return a matrix object:

    @classmethod
    def eval(cls, *_args):
        """."""

        if not _args:
            return

        u = _args[0]
        v = _args[1]

        uxv = Tuple(u[1]*v[2] - u[2]*v[1],
                     u[2]*v[0] - u[0]*v[2],
                     u[0]*v[1] - u[1]*v[0])
        return Matrix(uxv)

This is related to the new API design in SymPDE and in particular issue pyccel/sympde#150.

I did not check if this produces the correct results, i.e. if the transformation rule I applied is correct. But at least when comparing them on a domain without a mapping and with an identity mapping produces the same operator.

[jowezarek]

I wrote a test, trying to see if I can evaluate
$\int_{\Omega}\nabla\times(u\times B)\circ\nabla\times(v\times B)$
for
$B(x, y, z) = (0, 0, 1)^T,$
$u(x, y, z) = (-sin(x), -cos(z), 0)^T,$
$v(x, y, z) = (cos(x), sin(z), 0)^T.$
It then holds
$\nabla\times(u\times B)\circ\nabla\times(v\times B) = sin(x)cos(x) + sin(z)cos(z)$
with antiderivative
$\frac{1}{2}sin^2(x) + \frac{1}{2}sin^2(z).$

The corresponding integral over $\Omega=(0,1)^3$ should evaluate to $sin^2(1) \approx 0.7080734182735712$, e.g., in the case of no mapping, an identity mapping or a Collela mapping.

The corresponding integral over $\Omega=(0,2)^3$ should evaluate to $4sin^2(2) \approx 3.307287241727224$, e.g., in the case of a linear mapping of the form $(2x_1, 2x_2, 2x_3)$.

import  time
from    mpi4py   import MPI
import  numpy    as np

from    sympde.calculus                     import dot, inner, curl, cross
from    sympde.expr                         import BilinearForm, integral
from    sympde.topology                     import element_of, elements_of, Derham, Cube, Mapping
from    sympde.utilities.utils              import plot_domain

from    psydac.api.discretization           import discretize
from    psydac.api.settings                 import PSYDAC_BACKEND_GPYCCEL
from    psydac.linalg.solvers               import inverse

class IdentityMapping(Mapping):

    _expressions = {'x': 'x1',
                    'y': 'x2',
                    'z': 'x3'}
    
    _ldim        = 3
    _pdim        = 3

class LinearMapping(Mapping):

    _expressions = {'x': '2*x1',
                    'y': '2*x2',
                    'z': '2*x3'}
    
    _ldim        = 3
    _pdim        = 3

# valid??
class CollelaMapping3D(Mapping):
    _expressions = {'x': 'x1 + (a/2)*sin(2.*pi*(x1-0.5))*sin(2.*pi*(x2-0.5))',
                    'y': 'x2 + (a/2)*sin(2.*pi*(x1-0.5))*sin(2.*pi*(x2-0.5))',
                    'z': 'x3'}

    _ldim        = 3
    _pdim        = 3

mapping_options = [None, IdentityMapping('I'), LinearMapping('L'), CollelaMapping3D('C', a=0.3)]
# mapping_options = [CollelaMapping3D('C', a=0.3), ]

comm        = MPI.COMM_WORLD
mpi_rank    = comm.rank
backend     = PSYDAC_BACKEND_GPYCCEL

ncells      = [32, 32, 32] # [16, 16, 16]
degree      = [3, 3, 3] # [2, 2, 2]
periodic    = [False, False, False]

logical_domain = Cube('C', bounds1=(0, 1), bounds2=(0, 1), bounds3=(0, 1))

for i, mapping in enumerate(mapping_options):

    if mapping is not None:
        domain  = mapping(logical_domain)
    else:
        domain  = logical_domain

    #plot_domain(domain, draw=True, isolines=True)

    derham = Derham(domain)

    domain_h = discretize(domain, ncells=ncells, periodic=periodic, comm=comm)
    derham_h = discretize(derham, domain_h, degree=degree)

    # -----------------------------------------------------------------------------------------------------

    V1 = derham.V1
    V2 = derham.V2

    V1h = derham_h.V1
    V2h = derham_h.V2

    V1cs  = V1h.coeff_space
    V2cs  = V2h.coeff_space

    u1, v1, h   = elements_of(V1, names='u1, v1, h')
    B           = element_of (V2, name='B')

    # Hcurl mass matrix
    a1      = BilinearForm((u1, v1), integral(domain, dot(u1, v1)))
    # "q" for "historic" reasons, "n" as in "assemble every time step n"
    qn      = BilinearForm((u1, v1), integral(domain, inner( u1,           cross(B, v1))))
    # CuCu (curl-curl)
    cucu    = BilinearForm((u1, v1), integral(domain, inner(curl(u1), curl(v1))))
    # If successfully assembled: The entire preconditioner (pc) matrix (m)
    pcm     = BilinearForm((u1, v1), integral(domain, inner(curl(cross(u1, B)), curl(cross(v1, B)))))

    # We compare pcm with a combination of a1, qn and cucu



    t0 = time.time()
    a1_h    = discretize(a1, domain_h, (V1h, V1h), backend=backend)
    t1 = time.time()
    if mpi_rank == 0:
        print(f'Discretizing a1         : {(t1-t0):.3g} sec.')

    t0 = time.time()
    qn_h    = discretize(qn, domain_h, (V1h, V1h), backend=backend)
    t1 = time.time()
    if mpi_rank == 0:
        print(f'Discretizing qn         : {round(t1-t0, 2)} sec.')

    t0 = time.time()
    cucu_h  = discretize(cucu, domain_h, (V1h, V1h), backend=backend)
    t1 = time.time()
    if mpi_rank == 0:
        print(f'Discretizing cucu       : {round(t1-t0, 2)} sec.')

    t0 = time.time()
    pcm_h  = discretize(pcm, domain_h, (V1h, V1h), backend=backend)
    t1 = time.time()
    if mpi_rank == 0:
        print(f'Discretizing pca        : {round(t1-t0, 2)} sec.')

    M1   = a1_h.assemble()
    CUCU = cucu_h.assemble()

    _, P1, P2, _  = derham_h.projectors()

    G, C, D = derham_h.derivatives_as_matrices

    # We choose u, v, B such that
    # inner( curl(cross(u, B)), curl(cross(v, B)) )
    # = sin(z)cos(z) + sin(x)cos(x)

    # Integrated over (0,1)^3, this evaluates to sin(1)^2 = 0.7080734182735712
    # (no mapping, IdentityMapping, CollelaMapping)
    # Integrated over (0,2)^3, this evaluates to 4*sin(2)^2 = 3.307287241727224

    B1fun = lambda x, y, z: 0
    B2fun = lambda x, y, z: 0
    B3fun = lambda x, y, z: 1
    Bfun = (B1fun, B2fun, B3fun)

    u1fun = lambda x, y, z: -np.sin(x)
    u2fun = lambda x, y, z: -np.cos(z)
    u3fun = lambda x, y, z: 0
    ufun = (u1fun, u2fun, u3fun)

    v1fun = lambda x, y, z: np.cos(x)
    v2fun = lambda x, y, z: np.sin(z)
    v3fun = lambda x, y, z: 0
    vfun = (v1fun, v2fun, v3fun)

    BField = P2(Bfun)
    ucoeff = P1(ufun).coeffs
    vcoeff = P1(vfun).coeffs

    # 1.: M1_inv @ Qn is a L2 projection of cross(u, B) into H(curl)
    #     Combined with CUCU, this should mimic pcm.
    Qn      = qn_h.assemble(B=BField)
    M1_inv  = inverse(M1, 'cg', tol=1e-8, maxiter=1000)
    PC      = Qn.T @ M1_inv @ CUCU @ M1_inv @ Qn

    PCM     = pcm_h.assemble(B=BField)

    value   = PC.dot_inner(ucoeff, vcoeff)
    value2  = PCM.dot_inner(ucoeff, vcoeff)

    expected_value = 0.7080734182735712 if i != 2 else 3.307287241727224
    print(f'Mapping: {mapping} - Expected Value     : {expected_value}')
    print(f'PC L2 workaround value                  : {value}')
    print(f'PCM direct value                        : {value2}')
    print()

And while things worked out for no mapping and the identity mapping, the other two cases did not work out.

I changed two things in the code and now it seems like things are working. But I think we should do some more testing.

            elif isinstance(arg, Cross):

                arg = Cross(args[0].expr, args[1].expr)
                
                J   = mapping.jacobian

                #return J * curl( J.inv() * arg)
                return (J/J.det()) * curl( J.det() * J.inv() * arg )

I don't think there's a good reason for why we should be able to cancel out a 1/J.det() outside of the curl with a J.det() inside of the curl.
Then I had a look at ChatGPTs proposed pullback formula. At some point it wrote
$\widehat{u\times b} = det J J^{-1} (u\times b)$
which is the pullback to the Hdiv space, not the Hcurl space.
I think the correct way is
$\widehat{u\times b} = J^T (u\times b)$
and hence

            elif isinstance(arg, Cross):

                arg = Cross(args[0].expr, args[1].expr)
                
                J   = mapping.jacobian

                #return J * curl( J.inv() * arg)
                #return (J/J.det()) * curl( J.det() * J.inv() * arg )
                return (J/J.det()) * curl( J.T * arg )

Note that for complex mappings (e.g. Collela), discretizing the entire bilinear form
$(u, v)\mapsto\int_{\Omega}\nabla\times(u\times B)\circ\nabla\times(v\times B)$
requires a high resolution - in this example [32, 32, 32] Bspline cells and degree [3, 3, 3] Bsplines.
Choosing a worse resolution results in very wrong results. The workaround with the L2 projection performs better in that regard.