Transition from static to transient formulation

General
1

Hi DOLFINx community,

I am solving a transient EM axisymmetric magnetoquasistatic problem in DOLFINx using a backward Euler time discretization, and I am observing an unexpected behavior when transitioning from the static to the transient formulation.

Below I provide a minimal reproducible example consisting of:

  • a source coil driven by a prescribed current ramp (I0 + dIsdt * t), and

  • a passive conductive coil where eddy currents are induced.

By modifying dIsdt in a wide range (e.g. (0, 1e6]), one can observe that for sufficiently small ramp rates the eddy current evolution no longer follows the expected analytical RL curve.

More precisely:

  • For larger values of dIsdt, the transient solution is smooth and matches the analytical RL response very well.

  • As dIsdt decreases, a noticeable jump appears at the first transient time step.

  • After this first step, the solution evolves smoothly again, but the initial deviation remains.

This behavior seems to be directly linked to the first transient step (it is like the previous static result, prev_psi, not passed correctly to the transient analysis… indeed the issue is not present if I0 = 0).

Additional details:

  • DOLFINx version: 0.8.0

  • The problem is set up using LinearProblem

  • The system is fully reassembled at every time step (I am aware this is not optimal, but I intentionally started with a straightforward and “safe” implementation before optimizing the workflow)

I would greatly appreciate any suggestions to solve this problem.

Thanks in advance for your help.

import numpy as np
import gmsh
from mpi4py import MPI
from petsc4py import PETSc

import ufl
from ufl import TrialFunction, TestFunction, grad, dot, SpatialCoordinate

from dolfinx import fem
from dolfinx.fem import Function, functionspace, Constant, form
from dolfinx.fem.petsc import LinearProblem
from dolfinx.io import gmshio
from dolfinx.mesh import locate_entities_boundary


# -----------------------------------------------------------------------------
# 1. MESH GENERATION
# -----------------------------------------------------------------------------

gmsh.initialize()
gmsh.model.add("D_shape")

r_outer = 10.0
width = 0.1
shift = 1.0
mesh_size = 0.2

p1 = gmsh.model.occ.addPoint(0, r_outer, 0)
p2 = gmsh.model.occ.addPoint(0, 0, 0)
p3 = gmsh.model.occ.addPoint(0, -r_outer, 0)

arc = gmsh.model.occ.addCircleArc(p1, p2, p3)
line = gmsh.model.occ.addLine(p3, p1)

loop = gmsh.model.occ.addCurveLoop([arc, line])
d_surface = gmsh.model.occ.addPlaneSurface([loop])

square1 = gmsh.model.occ.addRectangle(shift, shift, 0, width, width)
square2 = gmsh.model.occ.addRectangle(shift, -shift, 0, width, width)

cut = gmsh.model.occ.cut([(2, d_surface)],
                         [(2, square1), (2, square2)],
                         removeObject=True,
                         removeTool=False)

gmsh.model.occ.synchronize()

gmsh.model.addPhysicalGroup(2, [cut[0][0][1]], 3)
gmsh.model.addPhysicalGroup(2, [square1], 1)
gmsh.model.addPhysicalGroup(2, [square2], 2)

gmsh.model.mesh.setSize(gmsh.model.getEntities(0), mesh_size)
gmsh.model.mesh.generate(2)

domain, cell_tags, facet_tags = gmshio.model_to_mesh(
    gmsh.model, MPI.COMM_WORLD, 0, gdim=2
)

gmsh.finalize()

dx = ufl.Measure("dx", domain=domain, subdomain_data=cell_tags)
tdim = domain.topology.dim

# -----------------------------------------------------------------------------
# 2. FUNCTION SPACES
# -----------------------------------------------------------------------------

V = functionspace(domain, ("Lagrange", 1))
Q = functionspace(domain, ("DG", 0))

psi = Function(V)
psi_prev = Function(V)
js = Function(Q)
sigma = Function(Q)

cells_square1 = cell_tags.find(1)
cells_square2 = cell_tags.find(2)

sigma.x.array[cells_square2] = 5.96e7
sigma.x.scatter_forward()

# -----------------------------------------------------------------------------
# 3. PHYSICAL PARAMETERS
# -----------------------------------------------------------------------------

MU0 = Constant(domain, PETSc.ScalarType(4e-7 * np.pi))
dt = Constant(domain, PETSc.ScalarType(1e-3))

I0 = 1e3
dIsdt = 1.0

R_self = 1.1077e-05
L_self = 3.9273e-06
tau_RL = L_self / R_self
dt.value = tau_RL / 10

# -----------------------------------------------------------------------------
# 4. BOUNDARY CONDITION
# -----------------------------------------------------------------------------

facets = locate_entities_boundary(
    domain, tdim - 1,
    lambda x: np.isclose(x[0], 0.0)
)

dofs = fem.locate_dofs_topological(V, tdim - 1, facets)
bc = fem.dirichletbc(PETSc.ScalarType(0.0), dofs, V)

# -----------------------------------------------------------------------------
# 5. STATIC INITIAL CONDITION
# -----------------------------------------------------------------------------

area_square1 = domain.comm.allreduce(
    fem.assemble_scalar(form(1.0 * dx(1))),
    op=MPI.SUM,
)

js.x.array[cells_square1] = I0 / area_square1
js.x.scatter_forward()

u = TrialFunction(V)
v = TestFunction(V)

x = SpatialCoordinate(domain)
r = x[0]
r_factor = 1 / (2.0 * np.pi * r)

a_static = (r_factor / MU0) * dot(grad(u), grad(v)) * dx
L_static = js * v * dx(1)

problem_static = LinearProblem(a_static, L_static, u=psi, bcs=[bc])
problem_static.solve()

psi_prev.x.array[:] = psi.x.array[:]
psi_prev.x.scatter_forward()
# -----------------------------------------------------------------------------
# 6. TIME LOOP (LinearProblem style)
# -----------------------------------------------------------------------------

t_end = 4 * tau_RL
t_values = np.arange(0, t_end + dt.value, dt.value)
Ieddy = np.zeros(len(t_values))

for i in range(1, len(t_values)):
    # Update source current
    js.x.array[cells_square1] = (
        (I0 + dIsdt * t_values[i]) / area_square1
    )
    js.x.scatter_forward()

    u = TrialFunction(V)
    v = TestFunction(V)

    # Base operators
    a = (r_factor / MU0) * dot(grad(u), grad(v)) * dx
    L = js * v * dx

    # Transient contribution (Backward Euler)
    a += (r_factor * sigma / dt) * u * v * dx
    L += (r_factor * sigma / dt) * psi_prev * v * dx

    problem = LinearProblem(a, L, u=psi, bcs=[bc])
    problem.solve()

    # Eddy current
    Je = fem.Function(Q)
    Je_expr = fem.Expression(
        -sigma * r_factor * (psi - psi_prev) / dt,
        Q.element.interpolation_points(),
    )
    Je.interpolate(Je_expr)

    Ieddy[i] = domain.comm.allreduce(
        fem.assemble_scalar(form(Je * dx(2))),
        op=MPI.SUM,
    )

    psi_prev.x.array[:] = psi.x.array[:]
    psi_prev.x.scatter_forward()

# -----------------------------------------------------------------------------
# 7. FINAL PLOT (FEM + ANALYTICAL)
# -----------------------------------------------------------------------------

if domain.comm.rank == 0:
    import matplotlib.pyplot as plt

    L_mutual = 1.62943168e-07

    Ieddy_analytic = (
        -L_mutual * dIsdt / R_self
        * (1 - np.exp(-t_values / tau_RL))
    )

    plt.figure(figsize=(8, 5))
    plt.plot(t_values, Ieddy, linewidth=2, label="DOLFINx FEM")
    plt.plot(t_values, Ieddy_analytic, "--", linewidth=2,
             label="Analytical RL")
    plt.xlabel("Time [s]")
    plt.ylabel("Eddy current [A]")
    plt.title("Transient Eddy Current Response")
    plt.legend()
    plt.grid(True)
    plt.tight_layout()
    plt.show()

Figure_1
Figure_1800×500 30 KB

Figure_2
Figure_2800×500 32.5 KB

2

I would additionally pass PETSc options to the solver to ensure that you are solving the problem to the expected precision:

petsc_options = {"ksp_monitor": None, "ksp_error_if_not_converged":True, "ksp_type":"preonly", "pc_type":"lu",
                 "pc_factor_mat_solver_type": "mumps"}
problem_static = LinearProblem(a_static, L_static, u=psi, bcs=[bc], petsc_options=petsc_options)

and similarly for the time dependent problem.

This does not resolve your problem for small values of dIsdt. However, I would like to the see the strong from of your equations, including initial conditions and boundary conditions.

To me it seems like your initial conditions (psi_prev=solution_of_static_psi) is not correct for the time-dependent problem, which would explain the jump at the first time step.

3

Thanks @dokken for the suggestion regarding the PETSc options. I added the solver configuration you proposed but, as you already stated, that doesn’t solve the jump in current.

The formulation is axisymmetric magnetoquasistatic, written in terms of the magnetic flux function \psi = 2 \pi r A_\phi. The strong form reads:

\frac{1}{\mu_0}\,\nabla \cdot \left( \frac{1}{2\pi r}\,\nabla \psi \right) + \sigma \,\frac{1}{2\pi r}\,\frac{\partial \psi}{\partial t} = j_s \tag{1}

in the computational domain.

  • In the source coil region: j_s \neq 0
  • In the passive conductive coil: \sigma > 0, j_s = 0
  • In the air region: \sigma = 0, j_s = 0

As boundary conditions:

  • Dirichlet condition \psi = 0 on r = 0
  • Dirichlet condition \psi = 0 on \partial \Omega
    The external boundary is placed sufficiently far from the sources so that the magnetic field can reasonably be assumed to vanish there. This is the same approximation commonly used also in other FEM codes.

The initial condition is obtained from a static solve

\frac{1}{\mu_0}\,\nabla \cdot \left( \frac{1}{2\pi r}\,\nabla \psi \right) = j_s(t=0) \tag{2}

so that

\psi^{0} = \psi_{\text{static}} \tag{3}

This is used as the initial state for the backward Euler scheme

\frac{1}{\mu_0}\,\nabla \cdot \left( \frac{1}{2\pi r}\,\nabla \psi^{n+1} \right) + \sigma \,\frac{1}{2\pi r}\, \frac{\psi^{n+1} - \psi^{n}}{\Delta t} = j_s^{n+1} \tag{4}

I do not see an obvious reason why the static initialisation should be inappropriate, but I may be overlooking something.

Thanks for your support.

4

I just noticed that, in the example, I forgot to impose the null boundary condition on the external boundary. Since the condition is also enforced on the axis, the following snippet should replace the boundary condition section in the example:

facets = locate_entities_boundary(
domain,
tdim - 1,
lambda x: np.full(x.shape[1], True)
)

dofs = fem.locate_dofs_topological(V, tdim - 1, facets)

bc = fem.dirichletbc(PETSc.ScalarType(0.0), dofs, V)

In any case, the results remain qualitatively unchanged.

5

I tried looking at this for a bit, but struggle to pinpoint the error. What I find odd, is that if you uncomment the update of js in the timeloop, then one expects a stagnant solution (IC satisfies steady PDE and no additional forcing), but you don’t. That’s a red flag to me, and that’s where I would continue the debugging.

Additional notes: a lot of the code inside of your time loop should be moved out. See here a reduced code with petsc options per @dokken’s suggestion:

import numpy as np
import gmsh
from mpi4py import MPI
from petsc4py import PETSc

import ufl
from ufl import TrialFunction, TestFunction, grad, dot, SpatialCoordinate

from dolfinx import fem
from dolfinx.fem import Function, functionspace, Constant, form
from dolfinx.fem.petsc import LinearProblem
from dolfinx.io import gmshio
from dolfinx.mesh import locate_entities_boundary


# -----------------------------------------------------------------------------
# 1. MESH GENERATION
# -----------------------------------------------------------------------------

gmsh.initialize()
gmsh.model.add("D_shape")

r_outer = 10.0
width = 0.1
shift = 1.0
mesh_size = 0.2

p1 = gmsh.model.occ.addPoint(0, r_outer, 0)
p2 = gmsh.model.occ.addPoint(0, 0, 0)
p3 = gmsh.model.occ.addPoint(0, -r_outer, 0)

arc = gmsh.model.occ.addCircleArc(p1, p2, p3)
line = gmsh.model.occ.addLine(p3, p1)

loop = gmsh.model.occ.addCurveLoop([arc, line])
d_surface = gmsh.model.occ.addPlaneSurface([loop])

square1 = gmsh.model.occ.addRectangle(shift, shift, 0, width, width)
square2 = gmsh.model.occ.addRectangle(shift, -shift, 0, width, width)

cut = gmsh.model.occ.cut([(2, d_surface)],
                         [(2, square1), (2, square2)],
                         removeObject=True,
                         removeTool=False)

gmsh.model.occ.synchronize()

gmsh.model.addPhysicalGroup(2, [cut[0][0][1]], 3)
gmsh.model.addPhysicalGroup(2, [square1], 1)
gmsh.model.addPhysicalGroup(2, [square2], 2)

gmsh.model.mesh.setSize(gmsh.model.getEntities(0), mesh_size)
gmsh.model.mesh.generate(2)

domain, cell_tags, facet_tags = gmshio.model_to_mesh(
    gmsh.model, MPI.COMM_WORLD, 0, gdim=2
)

gmsh.finalize()

dx = ufl.Measure("dx", domain=domain, subdomain_data=cell_tags)
tdim = domain.topology.dim

# -----------------------------------------------------------------------------
# 2. FUNCTION SPACES
# -----------------------------------------------------------------------------

V = functionspace(domain, ("Lagrange", 1))
Q = functionspace(domain, ("DG", 0))

psi = Function(V)
psi_prev = Function(V)
js = Function(Q)
sigma = Function(Q)

cells_square1 = cell_tags.find(1)
cells_square2 = cell_tags.find(2)

sigma.x.array[cells_square2] = 5.96e7
sigma.x.scatter_forward()

# -----------------------------------------------------------------------------
# 3. PHYSICAL PARAMETERS
# -----------------------------------------------------------------------------

MU0 = Constant(domain, PETSc.ScalarType(4e-7 * np.pi))
dt = Constant(domain, PETSc.ScalarType(1e-3))

I0 = 1e3
dIsdt = 1.0

R_self = 1.1077e-05
L_self = 3.9273e-06
tau_RL = L_self / R_self
dt.value = tau_RL / 10

# -----------------------------------------------------------------------------
# 4. BOUNDARY CONDITION
# -----------------------------------------------------------------------------

facets = locate_entities_boundary(
    domain, tdim - 1,
    lambda x: np.isclose(x[0], 0.0)
)

dofs = fem.locate_dofs_topological(V, tdim - 1, facets)
bc = fem.dirichletbc(PETSc.ScalarType(0.0), dofs, V)

# -----------------------------------------------------------------------------
# 5. STATIC INITIAL CONDITION
# -----------------------------------------------------------------------------

area_square1 = domain.comm.allreduce(
    fem.assemble_scalar(form(1.0 * dx(1))),
    op=MPI.SUM,
)

js.x.array[cells_square1] = I0 / area_square1
# js.x.scatter_forward() # I don't think you need this.

u = TrialFunction(V)
v = TestFunction(V)

x = SpatialCoordinate(domain)
r = x[0]
r_factor = 1 / (2.0 * np.pi * r)

a_static = (r_factor / MU0) * dot(grad(u), grad(v)) * dx
L_static = js * v * dx(1)

problem_static = LinearProblem(a_static, L_static, u=psi, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "mumps"})
problem_static.solve()

psi_prev.x.array[:] = psi.x.array[:]
psi_prev.x.scatter_forward()

# Transient contribution (Backward Euler)
a = a_static + (r_factor * sigma / dt) * u * v * dx
L = L_static + (r_factor * sigma / dt) * psi_prev * v * dx

problem = LinearProblem(a, L, u=psi, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "mumps"})

Je = fem.Function(Q)
Je_expr = fem.Expression(
    -sigma * r_factor * (psi - psi_prev) / dt,
    Q.element.interpolation_points(),
)

# -----------------------------------------------------------------------------
# 6. TIME LOOP (LinearProblem style)
# -----------------------------------------------------------------------------

t_end = 4 * tau_RL
t_values = np.arange(0, t_end + dt.value, dt.value)
Ieddy = np.zeros(len(t_values))

for i in range(1, len(t_values)):
    # Update source current
    js.x.array[cells_square1] = (
        (I0 + dIsdt * t_values[i]) / area_square1
    )
    # js.x.scatter_forward() # I don't think you need this

    problem.solve()

    # Eddy current
    Je.interpolate(Je_expr)

    Ieddy[i] = domain.comm.allreduce(
        fem.assemble_scalar(form(Je * dx(2))),
        op=MPI.SUM,
    )

    psi_prev.x.array[:] = psi.x.array[:]
    # psi_prev.x.scatter_forward() # I dont think you need this

# -----------------------------------------------------------------------------
# 7. FINAL PLOT (FEM + ANALYTICAL)
# -----------------------------------------------------------------------------

if domain.comm.rank == 0:
    import matplotlib.pyplot as plt

    L_mutual = 1.62943168e-07

    Ieddy_analytic = (
        -L_mutual * dIsdt / R_self
        * (1 - np.exp(-t_values / tau_RL))
    )

    plt.figure(figsize=(8, 5))
    plt.plot(t_values, Ieddy, linewidth=2, label="DOLFINx FEM")
    plt.plot(t_values, Ieddy_analytic, "--", linewidth=2,
             label="Analytical RL")
    plt.xlabel("Time [s]")
    plt.ylabel("Eddy current [A]")
    plt.title("Transient Eddy Current Response")
    plt.legend()
    plt.grid(True)
    plt.tight_layout()
    plt.show()

It looks to me like there are some ill-conditioning issues at play, with the extreme system parameters. In particular, the sigma.x.array[cells_square2] = 5.96e7 appears problematic, and sigma.x.array[:] = 5.96e7 already behaves a lot less unexpected.

I can’t really be of more help, I’m afraid.