Add PML demo - GSoC 2022

[mikics]

Currently, there are no demos in DOLFINx showing how to implement perfectly matched layers (PMLs) for time-harmonic electromagnetic problems.

The aim of this pull request is to add such demo. In particular, the demo will show how to implement rectangular PML for a simple electromagnetic scattering problem in two dimensions (GSoC 2022).

[jhale]

@mikics I'm getting a fail on this for complex64 in the efficiency assertions. In the eye norm the numerical calculations don't look too bad. Is the tolerance a bit tight?

[mikics]

@mikics I'm getting a fail on this for complex64 in the efficiency assertions. In the eye norm the numerical calculations don't look too bad. Is the tolerance a bit tight?

On my laptop (with complex32) the demo passes the test both in serial and in parallel (even though the error increases in the parallel version). I can try to refine the mesh around the boundary, and see if the error diminishes in the complex64 version.

[jhale]

Few things you could investigate:

• What's changed since 8d06b41 when you were running green?
• We're moving the images back to PETSc 3.17 this morning, had a lot of problems with 3.18. Worth waiting to see if that changes anything.
• 32 and 64 simply relate to the global indexing, not to the data type of the complex scalar value. I wouldn't expect to see any difference.
• The docker image fenicsproject/test-env:nightly has complex32 and complex64 versions of PETSc which you can build DOLFINx against and experiment yourself.

[jhale]

OK, it's running green now. Would still be interested to know if you can locate the cause of this error is before merging:

https://github.com/FEniCS/dolfinx/actions/runs/3733081716/jobs/6333395279

[mikics]

OK, it's running green now. Would still be interested to know if you can locate the cause of this error is before merging:

https://github.com/FEniCS/dolfinx/actions/runs/3733081716/jobs/6333395279

Will investigate that! We also have a duplicate module error now, since the same function for the calculation of analytical efficiencies is used for the demo for scattering boundary conditions. I will take care of that as well.

[jhale]

Thanks! Also, on this PR could you add doc links to the two that will be merged?

https://docs.fenicsproject.org/dolfinx/main/python/demos.html#all-demos

[mikics]

Few things you could investigate:

• What's changed since 8d06b41 when you were running green?

I have analyzed the ouputs of the CI checks for the 8d06b41 (working commit) and 0eca2f2 (first non-working commit). By looking at the petsc4py version, it seems that the two commits were checked with two different PETSc versions:

• 3.17.4 (working commit)
• 3.18.1 (not working commit)

I have also quickly checked the git diff of the two commits, and it doesn't seem there are major changes in the dolfinx code, but I have to better look into that to have a clear idea.

I would also like to check if the error obtained with PETSc 3.17.4 ("working" version) is affected in any way by the number of bits for the indexing, since the assert only tells us that the error is below 1%, but not if it changes. This may give us further hints on what is the problem.

[jhale]

Thanks for investigating. The problems related to PETSc 3.18 are more to do with how our CI system is setup (using different compilers for PETSc and DOLFINx leads to different macros being triggered related to complex types).

You are running green now on 3.17 so I think it's OK for us to merge.

[twmr]

Thx a lot for working on this! I would love to see a proper support for PMLs at some point in dolfinx. ngsolve e.g. has a SetPML method, which is documented here https://docu.ngsolve.org/latest/i-tutorials/unit-1.7-helmholtz/pml.html.

I have my own 2D version of an UPML, which I implemented in a slightly different way and therefore have some questions and feedback.

[twmr]

Here I would add the definition of x = SpatialCoordinate() and alpha and inline the pml_coordinates function. This IMO makes the code a bit easier to read, because you then don't have to scroll around in the notebook to understand how eps_*, mu_* are determined/defined.

This signure could be used:

def create_eps_mu(
        pml_x: bool,
        pml_y: bool,
        eps_bkg: Union[float, ufl.tensors.ListTensor],
        mu_bkg: Union[float, ufl.tensors.ListTensor],
        alpha: float=1,
) -> Tuple[ufl.tensors.ComponentTensor,
           ufl.tensors.ComponentTensor]:
    ...

eps_x, mu_x = create_eps_mu(True, False, eps_bkg, 1)
eps_y, mu_y = create_eps_mu(False, True, eps_bkg, 1)
eps_xy, mu_xy = create_eps_mu(True, True, eps_bkg, 1)

[mikics]

Do you mean something like this?

def create_eps_mu(
        domain,
        pml_x: bool,
        pml_y: bool,
        eps_bkg: Union[float, ufl.tensors.ListTensor],
        mu_bkg: Union[float, ufl.tensors.ListTensor],
        k0: float,
        l_dom: float,
        l_pml: float,
        alpha: float = 1) -> Tuple[ufl.tensors.ComponentTensor,
                                   ufl.tensors.ComponentTensor]:

    x = ufl.SpatialCoordinate(domain)

    if pml_x:
        x_pml = (x[0] + 1j * alpha / k0 * x
                 * (ufl.algebra.Abs(x[0]) - l_dom / 2)
                 / (l_pml / 2 - l_dom / 2)**2)
    else:
        x_pml = x[0]

    if pml_y:
        y_pml = (x[1] + 1j * alpha / k0 * x
                 * (ufl.algebra.Abs(x[1]) - l_dom / 2)
                 / (l_pml / 2 - l_dom / 2)**2)
    else:
        y_pml = x[1]

    pml_coordinates = ufl.as_vector((x_pml, y_pml))

    J = ufl.grad(pml_coordinates)

    J = ufl.as_matrix(((J[0, 0], 0, 0),
                       (0, J[1, 1], 0),
                       (0, 0, 1)))

    A = ufl.inv(J)
    eps_pml = ufl.det(J) * A * eps_bkg * ufl.transpose(A)
    mu_pml = ufl.det(J) * A * mu_bkg * ufl.transpose(A)

    return eps_pml, mu_pml

[twmr]

Yes exactly! I think the variable names (pml_x and pml_y) could be more descriptive, e.g., x_coord_in_pml and y_coord_in_pml, though.

[twmr]

Here you're defining three domains for the three parts of the PML, which you then use in your definition of the weak form. Was there a particular reason why you have done it like that instead of defining cellwise-constant fem.Functions over the entire domain? I'm asking because in my own implementation of an UPML, I'm using fem.Functions.

For eps you are even defining a fem.Function, but you are not using it in the PML region(s).

[mikics]

The main reason is that we wanted to explicitly show that we have three different transformations in the six different regions of the PML.

Besides, I think that using fem.Function would make the computation a bit slower, since you need to interpolate the function over a certain function space (but I might be wrong).

[twmr]

Besides, I think that using fem.Function would make the computation a bit slower, since you need to interpolate the function over a certain function space (but I might be wrong).

I agree with you that it would probably be slower, but I haven't verified if using fem.Function has a measurable impact. In my application the user has to provide a permittivity function and a permeability function for the entire domain (including the PML). The solver in my application is not aware that there is a PML in the geometry (In 1D the outgoing boundary condition can be imposed at the boundary of your geometry, i.e., a PML is not needed.). I think both approaches (yours and mine) have their pros and cons.

[twmr]

This is known as a UPML (uniaxial PML). Maybe this could be mentioned here somewhere.

[twmr]

For UPML implementations \varepsilon and \mu are given by

(http://www.mit.edu/~wsshin/pdf/shin2012jcp.pdf)

I think you don't have to use ufl.grad, but it is nice to see that this is possible using the ufl library.

[mikics]

Did not know about this definition of uniaxial PML, thank you!

I am not familiar with the difference between uniaxial PML and stretched-coordinate PML. The calculation reported in the demo are based on this work
, where the authors propose a complex coordinate transformation to implement PML for irregularly shaped domain.

In this demo the domain has a regular shape (a square), and therefore what you propose makes sense. However, with ufl.grad we are not constrained to any particular shape, and therefore I would say that this approach is "more general". For instance, in this other PR a spherical PML is expressed with the same formulation, and, in principle, it could be used for more complicated shapes as well.

[twmr]

Can't you use ufl.curl instead of defining your own 2D version of curl?

[mikics]

I think that using ufl.curl in this case is not possible since the problem is two-dimensional, and ufl.curl would try to evaluate the derivative of Es along the third dimension, resulting in an error. Therefore we need to define this version of the $\nabla \times$ operator.

[twmr]

At least for me it works in a 2D problem.

[mikics]

To be clearer, if I use this weak form:

F = - ufl.inner(ufl.curl(Es_3d), ufl.curl(v_3d)) * dDom \
    + eps * (k0**2) * ufl.inner(Es, v) * dDom \
    + (k0**2) * (eps - eps_bkg) * ufl.inner(Eb, v) * dDom \
    - ufl.inner(ufl.inv(mu_x) * ufl.curl(Es_3d), ufl.curl(v_3d)) * dPml_x \
    - ufl.inner(ufl.inv(mu_y) * ufl.curl(Es_3d), ufl.curl(v_3d)) * dPml_y \
    - ufl.inner(ufl.inv(mu_xy) * ufl.curl(Es_3d), ufl.curl(v_3d)) * dPml_xy \
    + (k0**2) * ufl.inner(eps_x * Es_3d, v_3d) * dPml_x \
    + (k0**2) * ufl.inner(eps_y * Es_3d, v_3d) * dPml_y \
    + (k0**2) * ufl.inner(eps_xy * Es_3d, v_3d) * dPml_xy

I get this error:

Traceback (most recent call last):
  File "/root/dolfinx/python/demo/demo_pml/demo_pml.py", line 444, in <module>
    a, L = ufl.lhs(F), ufl.rhs(F)
  File "/root/ufl/ufl/formoperators.py", line 69, in lhs
    form = expand_derivatives(form)
  File "/root/ufl/ufl/algorithms/ad.py", line 25, in expand_derivatives
    form = apply_algebra_lowering(form)
  File "/root/ufl/ufl/algorithms/apply_algebra_lowering.py", line 175, in apply_algebra_lowering
    return map_integrand_dags(LowerCompoundAlgebra(), expr)
  File "/root/ufl/ufl/algorithms/map_integrands.py", line 46, in map_integrand_dags
    return map_integrands(lambda expr: map_expr_dag(function, expr, compress),
  File "/root/ufl/ufl/algorithms/map_integrands.py", line 27, in map_integrands
    mapped_integrals = [map_integrands(function, itg, only_integral_type)
  File "/root/ufl/ufl/algorithms/map_integrands.py", line 27, in <listcomp>
    mapped_integrals = [map_integrands(function, itg, only_integral_type)
  File "/root/ufl/ufl/algorithms/map_integrands.py", line 35, in map_integrands
    return itg.reconstruct(function(itg.integrand()))
  File "/root/ufl/ufl/algorithms/map_integrands.py", line 46, in <lambda>
    return map_integrands(lambda expr: map_expr_dag(function, expr, compress),
  File "/root/ufl/ufl/corealg/map_dag.py", line 36, in map_expr_dag
    result, = map_expr_dags(function, [expression], compress=compress,
  File "/root/ufl/ufl/corealg/map_dag.py", line 99, in map_expr_dags
    r = handlers[v._ufl_typecode_](v, *[vcache[u] for u in v.ufl_operands])
  File "/root/ufl/ufl/algorithms/apply_algebra_lowering.py", line 168, in curl
    return as_vector((c(1, 2), c(2, 0), c(0, 1)))
  File "/root/ufl/ufl/algorithms/apply_algebra_lowering.py", line 161, in c
    return a[j].dx(i) - a[i].dx(j)
  File "/root/ufl/ufl/exproperators.py", line 498, in _dx
    d = Grad(d)
  File "/root/ufl/ufl/differentiation.py", line 140, in __new__
    dim = find_geometric_dimension(f)
  File "/root/ufl/ufl/domain.py", line 376, in find_geometric_dimension
    error("Cannot determine geometric dimension from expression.")
  File "/root/ufl/ufl/log.py", line 135, in error
    raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: Cannot determine geometric dimension from expression.

Maybe I just need to tweak something in the code, but I don't know what. If you can share your example with UPML maybe I can spot what I am missing.

[mikics]

Thx a lot for working on this! I would love to see a proper support for PMLs at some point in dolfinx. ngsolve e.g. has a SetPML method, which is documented here https://docu.ngsolve.org/latest/i-tutorials/unit-1.7-helmholtz/pml.html.

I have my own 2D version of an UPML, which I implemented in a slightly different way and therefore have some questions and feedback.

Thanks for giving precious comments on this PR! I will look at them right now.