N1curl interpolation issue for dipole sources

General
1

Hello,

I already posted before on this forum about interpolation issues with N1curl elements. I received advice that worked well for homogeneous fields, but I was not able to get the same precision with dipolar fields.

Using the same mesh, I have interpolated a homogeneous potential field and the associated magnetic field to a very satisfactory precision using both N1curl and Lagrange elements.

When interpolating homogeneous field potential in N1curl
	L2 norm of A_nedelec_DG - A_DG = 3.21e-17
	L2 norm of B_nedelec_DG - B_DG = 3.21e-17
When interpolating homogeneous field potential in Lagrange
	L2 norm of A_lagrange_DG - A_DG = 8.13e-19
	L2 norm of B_lagrange_DG - B_DG = 8.13e-19

However, the same method did not work for a dipolar field.

When interpolating dipole potential in N1curl
	 L2 norm of A_nedelec_DG - A_DG = 1.06e-02
	 L2 norm of B_nedelec_DG - B_DG = 1.06e-02
When interpolating dipole potential in Lagrange
	 L2 norm of A_lagrange_DG - A_DG = 1.22e-17
	 L2 norm of B_lagrange_DG - B_DG = 1.22e-17

Is there anything I am missing regarding N1curl interpolation ? Why is interpolation so bad for the dipole field/potential ?

I have put the errors I get with homogenous field potentials.

Here is a (fairly long) MWE exploring my interpolation process, for clarity :

A_DG : Potential interpolated from numpy, into DG
A_nedelec_DG : Potential interpolated from numpy, into N1curl, into DG
A_lagrange_DG : Potential interpolated from numpy, into Lagrange, into DG
B_DG : Field interpolated from numpy, into DG
B_nedelec_DG : Field interpolated from numpy, into N1curl, into DG
B_lagrange_DG : Field interpolated from numpy, into Lagrange, into DG
curl_A_nedelec_DG: curl of the potential, interpolated from numpy, into N1curl, into DG

"""MWE highlighting my issues with Nédélec interpolation."""
import gmsh
import numpy as np
import pandas as pd
import dolfinx
from dolfinx.fem import FunctionSpace, VectorFunctionSpace, Function
from mpi4py import MPI
import ufl

VACUUM_PERMEABILITY = 4 * np.pi * 1e-7
DIPOLE_MU = np.array([0, 0, 1])
DIPOLE_POSITION = np.array([0, 0, 0])


def generate_mesh():
    """Generate empty domain."""
    bkg_tag = 0
    sphere_tag = 1

    # Mesh definition
    gmsh.initialize()
    gmsh.option.setNumber("General.Terminal", 0)
    gmsh.model.add("MWE")

    gmsh.model.occ.addSphere(0, 0, 0, 1, tag=bkg_tag)
    gmsh.model.occ.addSphere(0, 0, 0, 1e-1, tag=sphere_tag)
    gmsh.model.occ.synchronize()

    gmsh.model.addPhysicalGroup(3, [1], tag=bkg_tag)
    gmsh.model.addPhysicalGroup(3, [2], tag=sphere_tag)

    gmsh.option.setNumber("Mesh.Algorithm", 2)
    gmsh.model.mesh.generate(3)
    gmsh.model.mesh.refine()
    gmsh.model.mesh.optimize("Netgen")
    # gmsh.fltk.run()

    mesh, _, _ = dolfinx.io.gmshio.model_to_mesh(gmsh.model, MPI.COMM_WORLD, 0, gdim=3)
    gmsh.finalize()

    return mesh



def numpy_potential():
    """Get numpy function for dipolar potential field."""

    def background_potential(x):
        """Return dipolar potential."""
        temp_mu_dipole = np.einsum("n,d->dn", np.ones(x.shape[1]), DIPOLE_MU)
        temp_position = np.einsum("n,d->dn", np.ones(x.shape[1]), DIPOLE_POSITION)

        r_dipole = x - temp_position
        norm_r_dipole = np.linalg.norm(r_dipole, axis=0)
        produit_vect = np.cross(temp_mu_dipole, r_dipole, axis=0)

        # Potential expression obtained from https://en.wikipedia.org/wiki/Magnetic_dipole
        bkg_potential = (
            VACUUM_PERMEABILITY * produit_vect / (4 * np.pi * norm_r_dipole**3)
        )

        return bkg_potential

    return background_potential


def numpy_to_nedelec(mesh, numpy_potential):
    """Nedelec expression from numpy."""
    # V = FunctionSpace(mesh, ("N1curl", 1))
    V = VectorFunctionSpace(mesh, ("Lagrange", 1))
    A_nedelec = Function(V)
    A_nedelec.interpolate(numpy_potential())
    A_nedelec.x.scatter_forward()
    return A_nedelec


def numpy_to_DG(mesh, numpy_potential):
    """DG expression from numpy."""
    V_DG = VectorFunctionSpace(mesh, ("DG", 1))
    A_DG = Function(V_DG)
    A_DG.interpolate(numpy_potential())
    A_DG.x.scatter_forward()
    return A_DG, V_DG


mesh = generate_mesh()
background_potential = numpy_potential()

A_nedelec = numpy_to_nedelec(mesh, numpy_potential)
A_DG, V_DG = numpy_to_DG(mesh, numpy_potential)
comm = V_DG.mesh.comm

A_error = Function(V_DG)

A_nedelec_DG = Function(V_DG)
A_nedelec_DG.interpolate(A_nedelec)

A_error.x.array[:] = A_nedelec_DG.x.array - A_DG.x.array

error_L2 = dolfinx.fem.form(ufl.inner(A_error, A_error) * ufl.dx)
print(
    f"L2 norm of A_nedelec_DG - A_DG = {np.real(np.sqrt(comm.allreduce(dolfinx.fem.assemble_scalar(error_L2), MPI.SUM))):.2e}"
)

####################################################################################################
####################################################################################################

curl_A_DG = Function(V_DG)
curl_A_DG.interpolate(
    dolfinx.fem.Expression(
        ufl.as_vector(
            (
                A_DG.sub(2).dx(1) - A_DG.sub(1).dx(2),  # Dz/Dy - Dy/Dz
                A_DG.sub(0).dx(2) - A_DG.sub(2).dx(0),  # Dx/Dz - Dz/Dx
                A_DG.sub(1).dx(0) - A_DG.sub(0).dx(1),  # Dy/Dx - Dx/Dy
            )
        ),
        V_DG.element.interpolation_points(),
    )
)

curl_A_nedelec_DG = Function(V_DG)
curl_A_nedelec_DG.interpolate(
    dolfinx.fem.Expression(
        ufl.as_vector(
            (
                A_nedelec_DG.sub(2).dx(1) - A_nedelec_DG.sub(1).dx(2),  # Dz/Dy - Dy/Dz
                A_nedelec_DG.sub(0).dx(2) - A_nedelec_DG.sub(2).dx(0),  # Dx/Dz - Dz/Dx
                A_nedelec_DG.sub(1).dx(0) - A_nedelec_DG.sub(0).dx(1),  # Dy/Dx - Dx/Dy
            )
        ),
        V_DG.element.interpolation_points(),
    )
)

curl_error = Function(V_DG)
curl_error.x.array[:] = curl_A_nedelec_DG.x.array - curl_A_DG.x.array
curl_error_L2 = dolfinx.fem.form(ufl.inner(curl_error, curl_error) * ufl.dx)
print(
    f"L2 norm of curl_A_nedelec_DG - curl_A_DG = {np.real(np.sqrt(comm.allreduce(dolfinx.fem.assemble_scalar(curl_error_L2), MPI.SUM))):.2e}"
)


def numpy_field():
    """Get numpy function for dipolar magnetic field."""

    def background_field(x):
        """Return dipolar field."""
        ######################################### DIPOLE #########################################
        temp_mu_dipole = np.einsum("n,d->dn", np.ones(x.shape[1]), DIPOLE_MU)
        temp_position = np.einsum("n,d->dn", np.ones(x.shape[1]), DIPOLE_POSITION)
        r_dipole = x - temp_position
        norm_r_dipole = np.linalg.norm(r_dipole, axis=0)
        dot_product = np.einsum("dn,dn->n", temp_mu_dipole, r_dipole)

        bkg_field = (VACUUM_PERMEABILITY / (4 * np.pi * norm_r_dipole**3)) * (
            (3 * r_dipole * dot_product / norm_r_dipole**2) - temp_mu_dipole
        )
        return bkg_field

    return background_field


background_field = numpy_field()

B_nedelec = numpy_to_nedelec(mesh, numpy_potential)
B_DG, _ = numpy_to_DG(mesh, numpy_potential)
comm = V_DG.mesh.comm

B_error = Function(V_DG)

B_nedelec_DG = Function(V_DG)
B_nedelec_DG.interpolate(B_nedelec)

B_error.x.array[:] = B_nedelec_DG.x.array - B_DG.x.array
error_L2 = dolfinx.fem.form(ufl.inner(B_error, B_error) * ufl.dx)
print(
    f"L2 norm of B_nedelec_DG - B_DG = {np.real(np.sqrt(comm.allreduce(dolfinx.fem.assemble_scalar(error_L2), MPI.SUM))):.2e}"
)
####################################################################################################
B_error.x.array[:] = B_nedelec_DG.x.array - curl_A_nedelec_DG.x.array
error_L2 = dolfinx.fem.form(ufl.inner(B_error, B_error) * ufl.dx)
print(
    f"L2 norm of B_nedelec_DG - curl_A_nedelec_DG = {np.real(np.sqrt(comm.allreduce(dolfinx.fem.assemble_scalar(error_L2), MPI.SUM))):.2e}"
)

B_error.x.array[:] = B_DG.x.array - curl_A_DG.x.array
error_L2 = dolfinx.fem.form(ufl.inner(B_error, B_error) * ufl.dx)
print(
    f"L2 norm of B_DG - curl_A_DG = {np.real(np.sqrt(comm.allreduce(dolfinx.fem.assemble_scalar(error_L2), MPI.SUM))):.2e}"
)


## When interpolating dipole potential in N1curl
# L2 norm of A_nedelec_DG - A_DG = 1.06e-02 <- why is this error so high 
# L2 norm of B_nedelec_DG - B_DG = 1.06e-02 <- why is this error so high 
# L2 norm of curl_A_nedelec_DG - curl_A_DG = 5.50e-01
# L2 norm of B_nedelec_DG - curl_A_nedelec_DG = 7.20e-03 
# L2 norm of B_DG - curl_A_DG = 5.54e-01

## When interpolating dipole potential in Lagrange
# L2 norm of A_lagrange_DG - A_DG = 1.22e-17
# L2 norm of B_lagrange_DG - B_DG = 1.22e-17
# L2 norm of curl_A_lagrange_DG - curl_A_DG = 6.12e-16
# L2 norm of B_lagrange_DG - curl_A_lagrange_DG = 5.54e-01
# L2 norm of B_DG - curl_A_DG = 5.54e-01

## When interpolating homogeneous field potential in N1curl
# L2 norm of A_nedelec_DG - A_DG = 3.21e-17
# L2 norm of B_nedelec_DG - B_DG = 3.21e-17
# L2 norm of curl_A_nedelec_DG - curl_A_DG = 1.22e-16
# L2 norm of B_nedelec_DG - curl_A_nedelec_DG = 2.01e-01 <- is curl introducing big errors ?
# L2 norm of B_DG - curl_A_DG = 2.01e-01 <- is curl introducing big errors ?

## When interpolating homogeneous field potential in Lagrange
# L2 norm of A_lagrange_DG - A_DG = 8.13e-19
# L2 norm of B_lagrange_DG - B_DG = 8.13e-19
# L2 norm of curl_A_lagrange_DG - curl_A_DG = 3.61e-17
# L2 norm of B_lagrange_DG - curl_A_lagrange_DG = 2.01e-01 <- is curl introducing big errors ?
# L2 norm of B_DG - curl_A_DG = 2.01e-01 <- is curl introducing big errors ?

Any help is appreciated ! Thank you :slight_smile:

2

Please try to reduce the example, as this is fairly lengthy and time comsuming to parse.

1 Like
3

Hello dokken, thank you for taking the time. You are right, sorry for the lenghty code, here is a shorter working example demonstrating my issue concerning the dipolar potential field :

"""MWE highlighting my issues with Nédélec interpolation."""
import gmsh
import numpy as np
import dolfinx
from dolfinx.fem import FunctionSpace, VectorFunctionSpace, Function
from mpi4py import MPI
import ufl

VACUUM_PERMEABILITY = 4 * np.pi * 1e-7
DIPOLE_MU = np.array([0, 0, 1])
DIPOLE_POSITION = np.array([0, 0, 0])


def generate_mesh():
    """Generate empty domain."""

    # Mesh definition
    gmsh.initialize()
    gmsh.option.setNumber("General.Terminal", 0)
    gmsh.model.add("MWE")

    gmsh.model.occ.addSphere(0, 0, 0, 1)
    gmsh.model.occ.addSphere(0, 0, 0, 1e-1)
    gmsh.model.occ.synchronize()

    gmsh.model.addPhysicalGroup(3, [1])
    gmsh.model.addPhysicalGroup(3, [2])

    gmsh.option.setNumber("Mesh.Algorithm", 2)
    gmsh.model.mesh.generate(3)
    gmsh.model.mesh.refine()
    gmsh.model.mesh.optimize("Netgen")

    mesh, _, _ = dolfinx.io.gmshio.model_to_mesh(gmsh.model, MPI.COMM_WORLD, 0, gdim=3)
    gmsh.finalize()

    return mesh


def numpy_potential():
    """Get numpy function for dipolar potential field."""

    def background_potential(x):
        """Return dipolar potential."""
        temp_mu_dipole = np.einsum("n,d->dn", np.ones(x.shape[1]), DIPOLE_MU)
        temp_position = np.einsum("n,d->dn", np.ones(x.shape[1]), DIPOLE_POSITION)

        r_dipole = x - temp_position
        norm_r_dipole = np.linalg.norm(r_dipole, axis=0)
        produit_vect = np.cross(temp_mu_dipole, r_dipole, axis=0)

        # Potential expression obtained from https://en.wikipedia.org/wiki/Magnetic_dipole
        bkg_potential = (
            VACUUM_PERMEABILITY * produit_vect / (4 * np.pi * norm_r_dipole**3)
        )

        return bkg_potential

    return background_potential


def numpy_to_nedelec(mesh, numpy_potential):
    """Nedelec expression from numpy."""
    V = FunctionSpace(mesh, ("N1curl", 1))
    # V = VectorFunctionSpace(mesh, ("Lagrange", 1))
    A_nedelec = Function(V)
    A_nedelec.interpolate(numpy_potential())
    A_nedelec.x.scatter_forward()
    return A_nedelec


def numpy_to_DG(mesh, numpy_potential):
    """DG expression from numpy."""
    V_DG = VectorFunctionSpace(mesh, ("DG", 1))
    A_DG = Function(V_DG)
    A_DG.interpolate(numpy_potential())
    A_DG.x.scatter_forward()
    return A_DG, V_DG

mesh = generate_mesh()

A_nedelec = numpy_to_nedelec(mesh, numpy_potential)
A_DG, V_DG = numpy_to_DG(mesh, numpy_potential)
A_nedelec_DG = Function(V_DG)
A_nedelec_DG.interpolate(A_nedelec)

comm = V_DG.mesh.comm

A_error = Function(V_DG)
A_error.x.array[:] = A_nedelec_DG.x.array - A_DG.x.array
error_L2 = dolfinx.fem.form(ufl.inner(A_error, A_error) * ufl.dx)
error_L2 = np.real(
    np.sqrt(comm.allreduce(dolfinx.fem.assemble_scalar(error_L2), MPI.SUM))
)
print(f"L2 norm of A_nedelec_DG - A_DG = {error_L2:.2e}")

## When interpolating dipole potential in N1curl
# L2 norm of A_nedelec_DG - A_DG = 1.06e-02 /!\

## When interpolating dipole potential in Lagrange
# L2 norm of A_lagrange_DG - A_DG = 1.22e-17

## When interpolating homogeneous field potential in N1curl
# L2 norm of A_nedelec_DG - A_DG = 3.21e-17

## When interpolating homogeneous field potential in Lagrange
# L2 norm of A_lagrange_DG - A_DG = 8.13e-19
4

What kind of mesh are you trying to create here? The inner and outer sphere is not connected at all.

The mesh is singular at the origin, and as you are not removing this, you are experiencing the effect of the singularity, and the fact that a DG-1 space and a first order Nedelec space uses different evaluation points.

1 Like
5

Thank you for the pointer. The spheres are indeed not connected. I have modified the mesh to be the simplest possible (as per MWE guidelines) :

def generate_mesh():
    """Generate empty domain."""

    # Mesh definition
    gmsh.initialize()
    gmsh.model.occ.addSphere(0, 0, 0, 1)
    gmsh.model.occ.synchronize()
    gmsh.model.addPhysicalGroup(3, [1])
    gmsh.option.setNumber("Mesh.Algorithm", 2)
    gmsh.model.mesh.generate(3)
    gmsh.model.mesh.optimize("Netgen")
    
    mesh, _, _ = dolfinx.io.gmshio.model_to_mesh(gmsh.model, MPI.COMM_WORLD, 0, gdim=3)
    gmsh.finalize()

    return mesh

and still experience poor precision. To avoid any singularity issue inside the mesh, I have also put the dipole source outside of the domain.

DIPOLE_POSITION = np.array([10, 10, 10])

I tried the same configuration with a cubic dolfinx mesh and had the same effect :

def generate_mesh():

    """Generate empty domain using dolfinx."""
    return dolfinx.mesh.create_unit_cube(
        MPI.COMM_WORLD, 5, 5, 5, dolfinx.mesh.CellType.tetrahedron
    )

Here are the exact metrics for reference, as you can see, precision between Lagrange and N1curl is consistent for homogeneous fields but not dipoles :

# GMSH 
# When interpolating dipole potential in N1curl
L2 norm of A_nedelec_DG - A_DG = 1.56e-12 /!\
# When interpolating dipole potential in Lagrange
L2 norm of A_lagrange_DG - A_DG = 7.16e-27
# When interpolating homogeneous field potential in N1curl
L2 norm of A_nedelec_DG - A_DG = 5.98e-15
# When interpolating homogeneous field potential in Lagrange
L2 norm of A_lagrange_DG - A_DG = 2.14e-16

# Dolfinx mesh
# When interpolating dipole potential in N1curl
L2 norm of A_nedelec_DG - A_DG = 6.90e-13 /!\
# When interpolating dipole potential in Lagrange
L2 norm of A_lagrange_DG - A_DG = 1.45e-27
# When interpolating homogeneous field potential in N1curl
L2 norm of A_nedelec_DG - A_DG = 1.28e-15
# When interpolating homogeneous field potential in Lagrange
L2 norm of A_lagrange_DG - A_DG = 8.49e-17

Also, I would assume that the difference in evaluation points for DG-1 and N1curl does not affect the interpolation process without the singularity you pointed out. Otherwise, how can I reduce this effect without designing a denser mesh ?

6

Where are the large errors here? If the L2 norm is 10^{-13} The integrated error is of order 10^{-25} (remove square root), which is way below machine precision.

1 Like
7

Thank you for the quick reply ! Yes, I rushed by setting the dipole too far away from the domain, making the norm of the incident potential really small. I have normalized the potential field so that the norm of A_{DG} is close to 1. Here is my updated and improved example, showcasing errors for N1curl that are significantly above machine precision. Lagrange interpolation error is negligible.

import gmsh
import numpy as np
import dolfinx
from dolfinx.fem import FunctionSpace, VectorFunctionSpace, Function
from mpi4py import MPI
import ufl

VACUUM_PERMEABILITY = 4 * np.pi * 1e-7
DIPOLE_MU = np.array([0, 0, 1]) * 3.8e7
DIPOLE_POSITION = np.array([1.5, 1.5, 1.5])

def generate_mesh():
    """Generate empty domain using dolfinx."""
    return dolfinx.mesh.create_unit_cube(
        MPI.COMM_WORLD, 5, 5, 5, dolfinx.mesh.CellType.tetrahedron
    )


def numpy_potential():
    """Get numpy function for dipolar potential field."""

    def background_potential(x):
        """Return dipolar potential."""
        temp_mu_dipole = np.einsum("n,d->dn", np.ones(x.shape[1]), DIPOLE_MU)
        temp_position = np.einsum("n,d->dn", np.ones(x.shape[1]), DIPOLE_POSITION)

        r_dipole = x - temp_position
        norm_r_dipole = np.linalg.norm(r_dipole, axis=0)
        produit_vect = np.cross(temp_mu_dipole, r_dipole, axis=0)

        # Potential expression obtained from https://en.wikipedia.org/wiki/Magnetic_dipole
        bkg_potential = (
            VACUUM_PERMEABILITY * produit_vect / (4 * np.pi * norm_r_dipole**3)
        )
        return bkg_potential

    return background_potential


def numpy_to_nedelec(mesh, numpy_potential):
    """Nedelec expression from numpy."""
    V = FunctionSpace(mesh, ("N1curl", 1))
    # V = VectorFunctionSpace(mesh, ("Lagrange", 1))
    A_nedelec = Function(V)
    A_nedelec.interpolate(numpy_potential())
    A_nedelec.x.scatter_forward()
    return A_nedelec


def numpy_to_DG(mesh, numpy_potential):
    """DG expression from numpy."""
    V_DG = VectorFunctionSpace(mesh, ("DG", 1))
    A_DG = Function(V_DG)
    A_DG.interpolate(numpy_potential())
    A_DG.x.scatter_forward()
    return A_DG, V_DG


mesh = generate_mesh()

A_nedelec = numpy_to_nedelec(mesh, numpy_potential)
A_DG, V_DG = numpy_to_DG(mesh, numpy_potential)
A_nedelec_DG = Function(V_DG)
A_nedelec_DG.interpolate(A_nedelec)

comm = V_DG.mesh.comm

A_error = Function(V_DG)
A_error.x.array[:] = A_nedelec_DG.x.array - A_DG.x.array
error_L2 = dolfinx.fem.form(ufl.inner(A_error, A_error) * ufl.dx)
error_L2 = np.real(comm.allreduce(dolfinx.fem.assemble_scalar(error_L2), MPI.SUM))

A_DG_norm =  dolfinx.fem.form(ufl.inner(A_DG, A_DG) * ufl.dx)
A_DG_norm = np.real(comm.allreduce(dolfinx.fem.assemble_scalar(A_DG_norm), MPI.SUM))

print(f"Squared L2 norm of A_DG = {A_DG_norm:.2e}")
print(f"Squared L2 norm of A_nedelec_DG - A_DG = {error_L2:.2e}")

The resulting errors :

# GMSH mesh
Squared  L2 norm of A_DG = 1.03e+00

# When interpolating dipole potential in N1curl
Squared L2 norm of A_nedelec_DG - A_DG = 1.29e-02 /!\

# When interpolating dipole potential in Lagrange
Squared L2 norm of A_lagrange_DG - A_DG = 1.57e-32

# Dolfinx mesh
Squared  L2 norm of A_DG = 1.35e+00

# When interpolating dipole potential in N1curl
Squared L2 norm of A_nedelec_DG - A_DG = 1.75e-02 /!\

# When interpolating dipole potential in Lagrange
Squared L2 norm of A_lagrange_DG - A_DG = 2.51e-32
8

I am not sure why you expect the two to be exactly the same.

You are interpolating a smooth function into a function space (the nedelec space) that is a sub-space of the DG space. It will of course not be able to approximate the function to the same degree as the DG function.

However, if you have a function in N1curl and interpolate it to DG, they should be exactly the same, ie.

import gmsh
import numpy as np
import dolfinx
from dolfinx.fem import FunctionSpace, VectorFunctionSpace, Function
from mpi4py import MPI
import ufl

VACUUM_PERMEABILITY = 4 * np.pi * 1e-7
DIPOLE_MU = np.array([0, 0, 1]) * 3.8e7
DIPOLE_POSITION = np.array([1.5, 1.5, 1.5])

def generate_mesh():
    """Generate empty domain using dolfinx."""
    return dolfinx.mesh.create_unit_cube(
        MPI.COMM_WORLD, 5, 5, 5, dolfinx.mesh.CellType.tetrahedron
    )


def numpy_potential():
    """Get numpy function for dipolar potential field."""

    def background_potential(x):
        """Return dipolar potential."""
        temp_mu_dipole = np.einsum("n,d->dn", np.ones(x.shape[1]), DIPOLE_MU)
        temp_position = np.einsum("n,d->dn", np.ones(x.shape[1]), DIPOLE_POSITION)

        r_dipole = x - temp_position
        norm_r_dipole = np.linalg.norm(r_dipole, axis=0)
        produit_vect = np.cross(temp_mu_dipole, r_dipole, axis=0)

        # Potential expression obtained from https://en.wikipedia.org/wiki/Magnetic_dipole
        bkg_potential = (
            VACUUM_PERMEABILITY * produit_vect / (4 * np.pi * norm_r_dipole**3)
        )
        return bkg_potential

    return background_potential


def numpy_to_nedelec(mesh, numpy_potential):
    """Nedelec expression from numpy."""
    V = FunctionSpace(mesh, ("N1curl", 1))
    # V = VectorFunctionSpace(mesh, ("Lagrange", 1))
    A_nedelec = Function(V)
    A_nedelec.interpolate(numpy_potential())
    A_nedelec.x.scatter_forward()
    return A_nedelec


def numpy_to_DG(mesh, numpy_potential):
    """DG expression from numpy."""
    V_DG = VectorFunctionSpace(mesh, ("DG", 1))
    A_DG = Function(V_DG)
    A_DG.interpolate(numpy_potential())
    A_DG.x.scatter_forward()
    return A_DG, V_DG


mesh = generate_mesh()

A_nedelec = numpy_to_nedelec(mesh, numpy_potential)
A_DG, V_DG = numpy_to_DG(mesh, numpy_potential)
A_nedelec_DG = Function(V_DG)
A_nedelec_DG.interpolate(A_nedelec)
A_DG.name = "DG"
A_nedelec_DG.name = "Nedelec"

comm = V_DG.mesh.comm

error_L2 = dolfinx.fem.form(ufl.inner(A_nedelec - A_nedelec_DG, A_nedelec-A_nedelec_DG) * ufl.dx)
error_L2 = np.real(comm.allreduce(dolfinx.fem.assemble_scalar(error_L2), MPI.SUM))


print(f"Squared L2 norm of A_nedelec_DG - A_Nedelec = {error_L2:.2e}")

return

Squared L2 norm of A_nedelec_DG - A_Nedelec = 6.19e-30
1 Like
9

Thank you for your time. Maybe I have a misunderstanding of N1curl space.
If the potential I am interpolating checks A_{potential} \in H(curl), which it is supposed to as it is a physical dipolar field, should it not be properly interpolated in N1curl ?

It would seem to me that any physical magnetic field, and their associated potential (for the Coulomb gauge), should be interpolated as precisely in DG space as in N1curl. In particular, if N1curl interpolation of a smooth H(curl) function leads to a loss of precision compared to Lagrange or DG, does it simply mean that solving electromagnetic problems in N1curl requires a very dense mesh and high orders of elements when compared to other functional spaces ?

Since my use case is a scattering problem, I want to take advantage of the properties of N1curl functions (interface, gauging, …), and avoid DG and Lagrange functional space for the solution.

How can I interpolate into or out of N1curl space to get my solution in a way that I can interpret externally without too many loss ?

Is it possible, for example, to export directly A_{nedelec} or curl(A_{nedelec}), into a numpy array ? Does it involve any DOF transformation to get nodal evaluation ?

10

If the potential is in N1Curl(\Omega) where \Omega is your discretized domain, it should be accurately interpolate it.

As the field dipole is surely continuous and a higher order function (given it is a function of the radius), you could simply represent it using ufl.SpatialCoordinate (and ufl.as_vector` to describe the dipole_mu and dipole_position) or use a higher order space to more accurately capture it.

A source on the RHS of your problem does not have to live in the same space as your unknown.
You could interpolate it into any space (that being higher order Lagrange, DG,N1Curl, 3rd order N1Curl) etc.

1 Like
11

Thank you for the suggestions. I have tried using ufl instead and the formulation is indeed simpler for the dipole, although I observe the same results. I guess the error I have been experimenting may not be related to interpolation but to the solving process instead ?

You are right that applying sources from a space that I can interpolate to without issue can be a good choice. However I still need my solution to live in N1curl(\Omega) for it to have the desired properties.

When I try to solve for either a homogeneous or dipolar field in vaccum, Lagrange outperforms N1curl, which can be expected, but N1curl has a high relative error that stays near constant with mesh density. This is consistent with my other simulatioIn empty space, this does not matter to me as the properties I’m after are useful mostly in simulations with multiple mediums.ns with these elements, but is not supposed to be the case. For MWE simplicity, I am applying the source in the same space as the solution. Modifying this does not result in significant changes.

Here are some graphs showing error evolution when comparing dolfinx mesh vs gmsh mesh, and Lagrange vs N1curl, for both dipolar and homogeneous fields.

relative_error_mesh_dipolar_field
relative_error_mesh_dipolar_field1920×968 69.4 KB

relative_error_mesh_homo_field
relative_error_mesh_homo_field1366×700 40.1 KB

Why is N1curl error so consistently high ?

Maybe this warrants a new thread, as my original issue was apparently misplaced. I am not sure if this is an an error with dolfinx, my implementation of it, or if it makes physical sense. No literature I found gave me any reason to believe that precision would not increase with mesh density with Nedelec elements.

Associated MWE comparing precision between Lagrange and N1curl for homogeneous field:

import dolfinx
from dolfinx.fem import FunctionSpace, VectorFunctionSpace, Function
from mpi4py import MPI
import ufl
from ufl import TestFunction, TrialFunction, curl, dx, inner, div
from dolfinx.fem.petsc import LinearProblem

APPLIED_FIELD = np.array([1, 2, 3])
N = 6


def generate_mesh_dolfinx():
    """Generate empty domain using dolfinx."""
    return dolfinx.mesh.create_unit_cube(
        MPI.COMM_WORLD, N, N, N, dolfinx.mesh.CellType.tetrahedron
    )


def on_boundary(x):
    return (
        np.isclose(x[0], 0)
        | np.isclose(x[1], 0)
        | np.isclose(x[2], 0)
        | np.isclose(x[0], 1)
        | np.isclose(x[1], 1)
        | np.isclose(x[2], 1)
    )


def numpy_potential():
    def background_potential(x):
        bkg_potential = np.array(
            [
                APPLIED_FIELD[1] * x[2] / 2 - APPLIED_FIELD[2] * x[1] / 2,
                APPLIED_FIELD[2] * x[0] / 2 - APPLIED_FIELD[0] * x[2] / 2,
                APPLIED_FIELD[0] * x[1] / 2 - APPLIED_FIELD[1] * x[0] / 2,
            ]
        )
        return bkg_potential
    return background_potential


def numpy_to_nedelec(mesh, numpy_potential):
    """Nedelec expression from numpy."""
    V_nedelec = FunctionSpace(mesh, ("N1curl", 1))
    A_nedelec = Function(V_nedelec)
    A_nedelec.interpolate(numpy_potential())
    A_nedelec.x.scatter_forward()
    return A_nedelec, V_nedelec


def numpy_to_lagrange(mesh, numpy_potential):
    """Lagrange expression from numpy."""
    V_lagrange = VectorFunctionSpace(mesh, ("Lagrange", 1))
    A_lagrange = Function(V_lagrange)
    A_lagrange.interpolate(numpy_potential())
    A_lagrange.x.scatter_forward()
    return A_lagrange, V_lagrange


def solve_problem(applied_function, function_space, mesh):
    u = TrialFunction(function_space)
    v = TestFunction(function_space)

    a = inner(curl(u), curl(v)) * dx
    a += inner(div(u), div(v)) * dx
    L = inner(curl(applied_function), curl(v)) * dx
    L += inner(div(applied_function), div(v)) * dx
    # Boundary Conditions
    mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim)
    boundary_facets = dolfinx.mesh.locate_entities_boundary(
        mesh, mesh.topology.dim - 1, on_boundary
    )
    boundary_dofs = dolfinx.fem.locate_dofs_topological(
        function_space, mesh.topology.dim - 1, boundary_facets
    )
    bc = dolfinx.fem.dirichletbc(applied_function, boundary_dofs)
    solution_function = Function(function_space)

    problem = LinearProblem(
        a,
        L,
        u=solution_function,
        bcs=[bc],
        petsc_options={
            "ksp_type": "gmres",
            "pc_type": "eisenstat",
            "ksp_atol": 1e-10,
            "ksp_max_it": 10000, #no max, as long a it hits convergence criterion
        },
    )
    problem.solve()
    return solution_function


def compute_norm(fem_function, comm):
    error = dolfinx.fem.form(ufl.inner(fem_function, fem_function) * ufl.dx)
    error = np.real(comm.allreduce(dolfinx.fem.assemble_scalar(error), MPI.SUM))
    return error


mesh = generate_mesh_dolfinx()
background_potential = numpy_potential()

A_nedelec, V_nedelec = numpy_to_nedelec(mesh, numpy_potential)
A_lagrange, V_lagrange = numpy_to_lagrange(mesh, numpy_potential)

A_s_nedelec = solve_problem(
    applied_function=A_nedelec, function_space=V_nedelec, mesh=mesh
)
A_s_lagrange = solve_problem(
    applied_function=A_lagrange, function_space=V_lagrange, mesh=mesh
)

# Error
comm = V_nedelec.mesh.comm

Asn_An = compute_norm(A_s_nedelec - A_nedelec, comm=comm)
Asl_Al = compute_norm(A_s_lagrange - A_lagrange, comm=comm)
Asl = compute_norm(A_s_lagrange, comm=comm)
Asn = compute_norm(A_s_nedelec, comm=comm)
Al = compute_norm(A_lagrange, comm=comm)
An = compute_norm(A_nedelec, comm=comm)

print(f"||A_s_nedelec||**2 = {Asn:.2e}")
print(f"||A_s_lagrange||**2 = {Asl:.2e}")

print(f"||A_lagrange||**2 = {Al:.2e}")
print(f"||A_nedelec||**2 = {An:.2e}")
# relative errors
print(f"||A_lagrange - A_s_lagrange||**2/||A_lagrange||**2 = {Asl_Al/Al:.2e}")
print(f"||A_nedelec - A_s_nedelec||**2/||A_nedelec||**2 = {Asn_An/An:.2e}")

Any help (litterature, tips, examples) is deeply appreciated !
Thank you !