Normalization of the weak formulation

mathematics
1

Hi All,

I am still working on solving the maxwell equations to analyse rotary electric machines.
When the machine saturates, the solver struggles to find the correct solutions. One plausible reason for this is that the residual norm is in the order of 10^-10 as is shown in a snippet of the logs below.


  0 KSP Residual norm 2.486670416291e-02 
  1 KSP Residual norm 1.380370820022e-02 
  2 KSP Residual norm 9.620016501752e-03 
  3 KSP Residual norm 6.778876375838e-03 
  4 KSP Residual norm 4.834307630619e-03 
  5 KSP Residual norm 3.854141821222e-03 
  6 KSP Residual norm 3.210644257159e-03 
  7 KSP Residual norm 2.709679289817e-03 
  8 KSP Residual norm 2.328649140280e-03 
  9 KSP Residual norm 1.876253299590e-03 
 10 KSP Residual norm 1.274454537374e-03 
 11 KSP Residual norm 8.325825129116e-04 
 12 KSP Residual norm 5.757867813002e-04 
 13 KSP Residual norm 4.327386318483e-04 
... 
849 KSP Residual norm 2.568153274638e-08 
850 KSP Residual norm 2.549669236354e-08 
851 KSP Residual norm 2.536519874296e-08 
852 KSP Residual norm 2.529548877872e-08 
853 KSP Residual norm 2.519325438527e-08 
854 KSP Residual norm 2.504843037186e-08 
855 KSP Residual norm 2.488862126196e-08 
856 KSP Residual norm 2.469410241949e-08 
Linear nls_solve_ solve converged due to CONVERGED_RTOL iterations 856 res rate 0.989229 R^2 0.900292
  Residual norms for nls_solve_ solve.
  0 KSP Residual norm 2.469410242046e-08 
  1 KSP Residual norm 2.462401061996e-08 
  2 KSP Residual norm 2.461975515437e-08 
  3 KSP Residual norm 2.461965025717e-08 
  4 KSP Residual norm 2.461553746786e-08 
  5 KSP Residual norm 2.460004576825e-08 
  6 KSP Residual norm 2.456632517194e-08 
  7 KSP Residual norm 2.453535644054e-08 
  8 KSP Residual norm 2.450601617182e-08 
  9 KSP Residual norm 2.445647418676e-08 
 10 KSP Residual norm 2.437671907680e-08 
 11 KSP Residual norm 2.428631000229e-08 
 12 KSP Residual norm 2.418247659123e-08 
 13 KSP Residual norm 2.407611880470e-08 
 14 KSP Residual norm 2.400868329143e-08 
 15 KSP Residual norm 2.397678288475e-08 
 16 KSP Residual norm 2.394553128623e-08 
 17 KSP Residual norm 2.382607443439e-08 
 18 KSP Residual norm 2.352944673248e-08 
 19 KSP Residual norm 2.283270699366e-08 
 20 KSP Residual norm 2.193691902575e-08 
 21 KSP Residual norm 2.051287366009e-08 
 22 KSP Residual norm 1.873364885805e-08 
 23 KSP Residual norm 1.695480400251e-08 
 24 KSP Residual norm 1.627732173885e-08 
 25 KSP Residual norm 1.606847683670e-08 
 26 KSP Residual norm 1.590508972070e-08 
 27 KSP Residual norm 1.577590815169e-08 
 28 KSP Residual norm 1.541737984068e-08 
 29 KSP Residual norm 1.501066093126e-08 
 30 KSP Residual norm 1.480147504859e-08 
 31 KSP Residual norm 1.466179262518e-08 
 32 KSP Residual norm 1.453793866409e-08 
 33 KSP Residual norm 1.439628591814e-08 
 34 KSP Residual norm 1.415070796630e-08 
 35 KSP Residual norm 1.391827913853e-08 
 36 KSP Residual norm 1.368773272826e-08 
 37 KSP Residual norm 1.344939502189e-08 
 38 KSP Residual norm 1.311544639084e-08 
 39 KSP Residual norm 1.275379848178e-08 
 40 KSP Residual norm 1.244343672985e-08 
 41 KSP Residual norm 1.219182473133e-08 
 42 KSP Residual norm 1.187614102958e-08 
 43 KSP Residual norm 1.144368391527e-08 
 44 KSP Residual norm 1.122369325172e-08 
 45 KSP Residual norm 1.114966097833e-08 
 46 KSP Residual norm 1.108994318831e-08 
 47 KSP Residual norm 1.106811955681e-08 
 48 KSP Residual norm 1.104313208692e-08 
 49 KSP Residual norm 1.102025641004e-08 
 50 KSP Residual norm 1.095850114207e-08 
 51 KSP Residual norm 1.089196583288e-08 
 52 KSP Residual norm 1.082586075693e-08 
 53 KSP Residual norm 1.078094078071e-08 
 54 KSP Residual norm 1.074012735346e-08 
 55 KSP Residual norm 1.070722893052e-08 
 56 KSP Residual norm 1.069240997411e-08 
 57 KSP Residual norm 1.068885775634e-08 
 58 KSP Residual norm 1.068607358146e-08 
 59 KSP Residual norm 1.068243642196e-08 
 60 KSP Residual norm 1.068148245188e-08 
 61 KSP Residual norm 1.068147964556e-08 
 62 KSP Residual norm 1.068142444109e-08 
 63 KSP Residual norm 1.068105006020e-08 
 64 KSP Residual norm 1.068054101022e-08 
 65 KSP Residual norm 1.067969562952e-08 
 66 KSP Residual norm 1.067750562117e-08 
 67 KSP Residual norm 1.066839825941e-08 
 68 KSP Residual norm 1.064962376047e-08 
 69 KSP Residual norm 1.063323753890e-08 
 70 KSP Residual norm 1.062285726263e-08 
 71 KSP Residual norm 1.061568506463e-08 
 72 KSP Residual norm 1.060217534543e-08 
 73 KSP Residual norm 1.058592620148e-08 
 74 KSP Residual norm 1.057428659266e-08 
 75 KSP Residual norm 1.054556246684e-08 
 76 KSP Residual norm 1.050702464670e-08 
 77 KSP Residual norm 1.043881850391e-08 
 78 KSP Residual norm 1.037598438684e-08 
 79 KSP Residual norm 1.032607036079e-08 
 80 KSP Residual norm 1.023824247547e-08 
 81 KSP Residual norm 1.009643034957e-08 
 82 KSP Residual norm 9.804912725869e-09 
Linear nls_solve_ solve converged due to CONVERGED_ATOL iterations 82 res rate 0.986886 R^2 0.866716
  Residual norms for nls_solve_ solve.
  0 KSP Residual norm 9.804912725714e-09 
Linear nls_solve_ solve converged due to CONVERGED_ATOL iterations 0 res rate 1. R^2 0.
  Residual norms for nls_solve_ solve.
  0 KSP Residual norm 9.804912725714e-09 
Linear nls_solve_ solve converged due to CONVERGED_ATOL iterations 0 res rate 1. R^2 0.
  Residual norms for nls_solve_ solve.
  0 KSP Residual norm 9.804912725714e-09 
Linear nls_solve_ solve converged due to CONVERGED_ATOL iterations 0 res rate 1. R^2 0.
  Residual norms for nls_solve_ solve.
  0 KSP Residual norm 9.804912725714e-09 
Linear nls_solve_ solve converged due to CONVERGED_ATOL iterations 0 res rate 1. R^2 0.
  Residual norms for nls_solve_ solve.
  0 KSP Residual norm 9.804912725714e-09 
Linear nls_solve_ solve converged due to CONVERGED_ATOL iterations 0 res rate 1. R^2 0.
  Residual norms for nls_solve_ solve.
  0 KSP Residual norm 9.804912725714e-09 
Linear nls_solve_ solve converged due to CONVERGED_ATOL iterations 0 res rate 1. R^2 0.
  Residual norms for nls_solve_ solve.
  0 KSP Residual norm 9.804912725714e-09 
Linear nls_solve_ solve converged due to CONVERGED_ATOL iterations 0 res rate 1. R^2 0.
  Residual norms for nls_solve_ solve.
  0 KSP Residual norm 9.804912725714e-09 
Linear nls_solve_ solve converged due to CONVERGED_ATOL iterations 0 res rate 1. R^2 0.
  Residual norms for nls_solve_ solve.
  0 KSP Residual norm 9.804912725714e-09 
Linear nls_solve_ solve converged due to CONVERGED_ATOL iterations 0 res rate 1. R^2 0.
  Residual norms for nls_solve_ solve.
  0 KSP Residual norm 9.804912725714e-09 
Linear nls_solve_ solve converged due to CONVERGED_ATOL iterations 0 res rate 1. R^2 0.
  Residual norms for nls_solve_ solve.
  0 KSP Residual norm 9.804912725714e-09 
Linear nls_solve_ solve converged due to CONVERGED_ATOL iterations 0 res rate 1. R^2 0.
  Residual norms for nls_solve_ solve.
  0 KSP Residual norm 9.804912725714e-09 
Linear nls_solve_ solve converged due to CONVERGED_ATOL iterations 0 res rate 1. R^2 0.
  Residual norms for nls_solve_ solve.
  0 KSP Residual norm 9.804912725714e-09 
Linear nls_solve_ solve converged due to CONVERGED_ATOL iterations 0 res rate 1. R^2 0.
  Residual norms for nls_solve_ solve.
  0 KSP Residual norm 9.804912725714e-09 
Linear nls_solve_ solve converged due to CONVERGED_ATOL iterations 0 res rate 1. R^2 0.
  Residual norms for nls_solve_ solve.
  0 KSP Residual norm 9.804912725714e-09 
Linear nls_solve_ solve converged due to CONVERGED_ATOL iterations 0 res rate 1. R^2 0.
  Residual norms for nls_solve_ solve.
  0 KSP Residual norm 9.804912725714e-09 
Linear nls_solve_ solve converged due to CONVERGED_ATOL iterations 0 res rate 1. R^2 0.
  Residual norms for nls_solve_ solve.
  0 KSP Residual norm 9.804912725714e-09 
Linear nls_solve_ solve converged due to CONVERGED_ATOL iterations 0 res rate 1. R^2 0.
  Residual norms for nls_solve_ solve.
  0 KSP Residual norm 9.804912725714e-09 
Linear nls_solve_ solve converged due to CONVERGED_ATOL iterations 0 res rate 1. R^2 0.
  Residual norms for nls_solve_ solve.
  0 KSP Residual norm 9.804912725714e-09 
Linear nls_solve_ solve converged due to CONVERGED_ATOL iterations 0 res rate 1. R^2 0.
  Residual norms for nls_solve_ solve.
  0 KSP Residual norm 9.804912725714e-09 
Linear nls_solve_ solve converged due to CONVERGED_ATOL iterations 0 res rate 1. R^2 0.
  Residual norms for nls_solve_ solve.
  0 KSP Residual norm 9.804912725714e-09 
Linear nls_solve_ solve converged due to CONVERGED_ATOL iterations 0 res rate 1. R^2 0.
  Residual norms for nls_solve_ solve.
  0 KSP Residual norm 9.804912725714e-09 
Linear nls_solve_ solve converged due to CONVERGED_ATOL iterations 0 res rate 1. R^2 0.
  Residual norms for nls_solve_ solve.
  0 KSP Residual norm 9.804912725714e-09 
Linear nls_solve_ solve converged due to CONVERGED_ATOL iterations 0 res rate 1. R^2 0.
  Residual norms for nls_solve_ solve.
  0 KSP Residual norm 9.804912725714e-09 
Linear nls_solve_ solve converged due to CONVERGED_ATOL iterations 0 res rate 1. R^2 0.
  Residual norms for nls_solve_ solve.
  0 KSP Residual norm 9.804912725714e-09 
Linear nls_solve_ solve converged due to CONVERGED_ATOL iterations 0 res rate 1. R^2 0.
  Residual norms for nls_solve_ solve.
  0 KSP Residual norm 9.804912725714e-09 
Linear nls_solve_ solve converged due to CONVERGED_ATOL iterations 0 res rate 1. R^2 0.
  Residual norms for nls_solve_ solve.
  0 KSP Residual norm 9.804912725714e-09 
Linear nls_solve_ solve converged due to CONVERGED_ATOL iterations 0 res rate 1. R^2 0.
  Residual norms for nls_solve_ solve.
  0 KSP Residual norm 9.804912725714e-09 
Linear nls_solve_ solve converged due to CONVERGED_ATOL iterations 0 res rate 1. R^2 0.
  Residual norms for nls_solve_ solve.
  0 KSP Residual norm 9.804912725714e-09 

2025-12-04 14:50:49 - MATE - WARNING - -- Failed to solve the nonlinear problem: Newton solver did not converge because maximum number of iterations reached --

To improve this, I heard that it can be useful to normalize the weak form, or even define a dimensionless version. I can however, find very little information on this in literature, specifically concerning the maxwell equations.

Is this indeed a practice that is often used to improve solver robustness? And are there straightforward ways with which this normalization is usually performed? Or is a dimensionless formulation the way to go?

2

Can you share an MWE that reproduces your problem? Or link to the earlier thread if you already shared it before?

I myself have faced convergence issues modeling nonlinear magnetic materials. Generally, using linear solution as initial guess for the nonlinear problem has worked for me. The strategy did not work in some extreme cases. I could not try it in fenics models but generally in these situations I resort to iteratively ramping up the load using the preceding nonlinear result as initial guess for the new problem.

3

Simula published a book on the topic.

2 Likes
4

Thank you for your reply @sbhasan .

I’ve tried using the linear solution as an initial guess and load stepping and while it certainly improved things a bit, it does not seem to be sufficient. Below are two scripts, one to create a dummy machine geometry (run this first to create the mesh) and then another to perform the actual analysis. I hope you can shed some insight into why the solver is not converging!

import gmsh
import sys
import numpy as np
import pickle

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

mm = 1e-3

# Radii
r_rotor = 20*mm
r_airgap_in = 20*mm
r_airgap_out = 21*mm
r_stator = 40*mm
r_outer_air = 60*mm

# Counts
p_pairs = 4
poles = 2 * p_pairs
magnets = 8
slots = 48

mesh_size = 1e-3

# Center point for arcs
center_pt = gmsh.model.occ.addPoint(0, 0, 0)

def add_sector(r1, r2, a1, a2):
    p = [
        gmsh.model.occ.addPoint(r1*np.cos(a1), r1*np.sin(a1), 0, mesh_size),
        gmsh.model.occ.addPoint(r2*np.cos(a1), r2*np.sin(a1), 0, mesh_size),
        gmsh.model.occ.addPoint(r2*np.cos(a2), r2*np.sin(a2), 0, mesh_size),
        gmsh.model.occ.addPoint(r1*np.cos(a2), r1*np.sin(a2), 0, mesh_size)
    ]
    l1 = gmsh.model.occ.addLine(p[0], p[1])
    l2 = gmsh.model.occ.addCircleArc(p[1], center_pt, p[2])
    l3 = gmsh.model.occ.addLine(p[2], p[3])
    l4 = gmsh.model.occ.addCircleArc(p[3], center_pt, p[0])
    cl = gmsh.model.occ.addCurveLoop([l1, l2, l3, l4])
    return gmsh.model.occ.addPlaneSurface([cl])

# --- Airgap points ---
r_in = r_airgap_in
r_out = r_airgap_out
n_points = 256

inner_pts = []
outer_pts = []
for i in range(n_points):
    theta = 2*np.pi*i/n_points
    inner_pts.append(gmsh.model.occ.addPoint(r_in*np.cos(theta), r_in*np.sin(theta), 0, mesh_size))
    outer_pts.append(gmsh.model.occ.addPoint(r_out*np.cos(theta), r_out*np.sin(theta), 0, mesh_size))

inner_arcs = []
outer_arcs = []
for i in range(n_points):
    next_i = (i+1) % n_points
    inner_arcs.append(gmsh.model.occ.addCircleArc(inner_pts[i], center_pt, inner_pts[next_i]))
    outer_arcs.append(gmsh.model.occ.addCircleArc(outer_pts[i], center_pt, outer_pts[next_i]))

cl_outer = gmsh.model.occ.addCurveLoop(outer_arcs)
cl_inner = gmsh.model.occ.addCurveLoop(inner_arcs)
airgap_surf = gmsh.model.occ.addPlaneSurface([cl_outer, cl_inner])

# --- Rotor using inner_pts (shared with airgap) ---
rotor_arcs = []
rotor_lines = []
for i in range(n_points):
    next_i = (i+1) % n_points
    rotor_arcs.append(gmsh.model.occ.addCircleArc(inner_pts[i], center_pt, inner_pts[next_i]))
cl_rotor = gmsh.model.occ.addCurveLoop(rotor_arcs)
rotor = gmsh.model.occ.addPlaneSurface([cl_rotor])

# --- Stator ring with shared airgap and refined outer edge ---
# Outer points for refined mesh
n_outer = 512
stator_outer_points = []
stator_outer_arcs = []
for i in range(n_outer):
    theta = 2 * np.pi * i / n_outer
    stator_outer_points.append(gmsh.model.occ.addPoint(r_stator*np.cos(theta), r_stator*np.sin(theta), 0, mesh_size))

for i in range(n_outer):
    next_i = (i+1) % n_outer
    stator_outer_arcs.append(gmsh.model.occ.addCircleArc(stator_outer_points[i], center_pt, stator_outer_points[next_i]))

# Inner arcs: reuse outer arcs of airgap
stator_inner_arcs = outer_arcs  # outer_arcs from airgap

# Create curve loops
cl_inner = gmsh.model.occ.addCurveLoop(stator_inner_arcs)
cl_outer = gmsh.model.occ.addCurveLoop(stator_outer_arcs)

# Create plane surface with inner hole
stator_surf = gmsh.model.occ.addPlaneSurface([cl_outer, cl_inner])


# --- Magnets (8) ---
mag_surfs = []
for k in range(magnets):
    a1 = 2*np.pi*k/magnets - np.pi/magnets*0.6
    a2 = 2*np.pi*k/magnets + np.pi/magnets*0.6
    mag_surfs.append(add_sector(r_rotor*0.5, r_rotor*0.95, a1, a2))

# --- Coils (48 slots) ---
coil_surfs = []
for k in range(slots):
    a1 = 2*np.pi*k/slots + 0.01
    a2 = 2*np.pi*(k+1)/slots - 0.01
    coil_surfs.append(add_sector(r_airgap_out*1.05, r_stator*0.95, a1, a2))

gmsh.model.occ.synchronize()

# --- Boolean cuts, keep tools ---
rotor_cut, _ = gmsh.model.occ.cut(
    [(2, rotor)], [(2, m) for m in mag_surfs],
    removeObject=True, removeTool=False
)
rotor = rotor_cut[0][1]

stator_cut, _ = gmsh.model.occ.cut(
    [(2, stator_surf)], [(2, c) for c in coil_surfs],
    removeObject=True, removeTool=False
)
stator_ids = [s[1] for s in stator_cut]

gmsh.model.occ.synchronize()

# --- Physical groups ---
gmsh.model.addPhysicalGroup(2, [rotor], 1)
gmsh.model.setPhysicalName(2, 1, "Rotor")

gmsh.model.addPhysicalGroup(2, [airgap_surf], 2)
gmsh.model.setPhysicalName(2, 2, "Airgap")

gmsh.model.addPhysicalGroup(2, stator_ids, 3)
gmsh.model.setPhysicalName(2, 3, "Stator")

# --- Coils: split into consecutive pairs ---
phase_names = ["A", "B", "C"]
phase_groups = {"A": [], "B": [], "C": []}

# Loop over all 48 coils
for i, cid in enumerate(coil_surfs):
    # Determine phase: A=0, B=1, C=2 repeating every 6 coils
    block = (i // 2) % 3
    phase = phase_names[block]
    phase_groups[phase].append(cid)

# Now create physical groups of two neighboring coils
for phase in phase_names:
    coils = phase_groups[phase]
    n_pairs = len(coils) // 2
    for j in range(n_pairs):
        group_name = f"Phase{phase}{j+1}"
        group_coils = coils[2*j:2*j+2]
        gmsh.model.addPhysicalGroup(2, group_coils, 200 + j + (phase_names.index(phase)*10))
        gmsh.model.setPhysicalName(2, 200 + j + (phase_names.index(phase)*10), group_name)


# Individual magnets
for i, mid in enumerate(mag_surfs):
    tag = 100 + i
    gmsh.model.addPhysicalGroup(2, [mid], tag)
    gmsh.model.setPhysicalName(2, tag, f"Magnet_{i+1}")

# --- Build geometry dictionary ---
geometry = {}

def geom_entry(name, mat):
    phys = gmsh.model.getPhysicalGroups()
    for dim, tag in phys:
        if gmsh.model.getPhysicalName(dim, tag) == name:
            entities = gmsh.model.getEntitiesForPhysicalGroup(dim, tag)
            return {"geom": {"name": name, "tag": tag, "surface": entities},
                    "material": mat}
    raise RuntimeError(f"Physical group {name} not found")

geometry["Center"] = {"center": 1}
geometry["Rotor"]  = geom_entry("Rotor",  {"mu_r": 2000})
geometry["Stator"] = geom_entry("Stator", {"mu_r": 2000})
geometry["Airgap"] = geom_entry("Airgap", {"mu_r": 1.0})
# --- Geometry dictionary for sub-phase coils ---
for phase in ["A", "B", "C"]:
    for i in range(8):
        name = f"Phase{phase}{i+1}"
        geometry[name] = geom_entry(name, {"mu_r": 1.0})

for i in range(magnets):
    name = f"Magnet_{i+1}"
    geometry[name] = geom_entry(name, {"mu_r": 1.05})

with open("geometry.pkl", "wb") as f:
    pickle.dump(geometry, f)

# --- Mesh ---
gmsh.model.mesh.generate(2)
gmsh.write("simple_ipm_machine.msh")

if "-nopopup" not in sys.argv:
    gmsh.fltk.run()

gmsh.finalize()

Analysis:

import pickle

import matplotlib
import numpy as np
import pyvista
import ufl
from basix.ufl import element
from dolfinx import default_scalar_type
from dolfinx.fem import functionspace, locate_dofs_topological, dirichletbc, Function, Expression, form, assemble_scalar
from dolfinx.io.gmshio import read_from_msh
from dolfinx.mesh import locate_entities_boundary
from dolfinx.plot import vtk_mesh
from mpi4py.MPI import COMM_WORLD
from petsc4py import PETSc
from scipy.interpolate import PchipInterpolator
from ufl import Measure, TrialFunction, TestFunction, as_ufl, inner, curl, grad

I_ph = 200*9
position = 0.0

bh_curve =  [(0.0000, 0.0000), (0.2447, 38.9522), (0.5008, 60.7621), (0.7509, 104.1758), (1.0000, 187.0970), (1.2522, 426.7090),
            (1.3949, 1157.2412), (1.4994, 2804.8711), (1.5934, 5220.9581), (1.7116, 10592.0), (1.7507, 13305.0),
            (1.7917, 16906.0), (1.8303, 21753.0), (1.8486, 25194.0), (1.8664, 29701.0), (1.8883, 37281.0),
            (1.9057, 45251.0), (1.9230, 55000.0), (1.9404, 66226.0), (1.9578, 78595.0), (1.9751, 91776.0),
            (2.0055, 115790.0), (2.1117, 200310.0), (2.4947, 505020.0), (3.5124, 1314900.0)]

def nu(u: Function) -> Function:
    from scipy.interpolate import CubicSpline
    mu0 = 4 * np.pi * 1e-7  # Permeability of free space
    B_vals, H_vals = zip(*bh_curve)
    B_vals = np.array([b for b in B_vals])
    H_vals = np.array([h for h in H_vals])
    y = np.diff(np.insert(H_vals, 0, -H_vals[1])) / np.diff(np.insert(B_vals, 0, -B_vals[1]))
    y = np.clip(y, a_min=1/7550 * 1/mu0, a_max=None)

    # cs = CubicSpline(B_vals, y, bc_type='clamped')
    cs = PchipInterpolator(B_vals, y)
    B_breaks, B_coefs = cs.x, cs.c

    dm = V_DG.dofmap
    c2d = dm.cell_dofs

    nu_f = Function(V_DG)
    nu_f.x.array[:] = 1 / mu0
    nu_func = Function(V_DG)
    nu_func.x.array[:] = 1 / mu0

    B_norm = ufl.sqrt(u.dx(1) ** 2 + (-u.dx(0)) ** 2)  / ufl.sqrt(2)

    nu_interp = 0
    for i in range(len(B_breaks) - 1):
        # mask = (B_norm >= B_breaks[i]) & (B_norm < B_breaks[i + 1])
        delta_x = B_norm - B_breaks[i]
        if i == len(B_breaks) - 2:
            # Completely saturated
            poly = 1 / mu0
            mask = B_norm >= B_breaks[i]
        else:
            poly = (B_coefs[0, i] * delta_x ** 3 +
                    B_coefs[1, i] * delta_x ** 2 +
                    B_coefs[2, i] * delta_x +
                    B_coefs[3, i])
            mask = ufl.And(B_norm >= B_breaks[i], B_norm < B_breaks[i + 1])
        nu_interp += ufl.conditional(mask, poly, 0)
    nu_expr = Expression(nu_interp, V_DG.element.interpolation_points())
    nu_func.interpolate(nu_expr)

    for key in geometry.keys() - {"Center"}:
        out = geometry[key]
        tag = out["geom"]["tag"]

        if key in ["Stator", "Rotor"]:
            for c in ct.find(tag):
                nu_f.x.array[c2d(c)] = nu_func.x.array[c2d(c)]
        else:
            mur = out["material"]["mu_r"]
            nu_val = 1 / (mu0 * mur)
            for c in ct.find(tag):
                nu_f.x.array[c2d(c)] = nu_val
    nu_f.x.scatter_forward()
    return nu_f

# Load the geometry pickle
with open("geometry.pkl", 'rb') as f:
    geometry = pickle.load(f)
# mesh, ct, ft = read_from_msh("full_geometry.msh", comm=COMM_WORLD, rank=0, gdim=2)
mesh, ct, ft = read_from_msh("simple_ipm_machine.msh", comm=COMM_WORLD, rank=0, gdim=2)
mesh.topology.create_connectivity(2, 1)
mesh.topology.create_connectivity(1, 2)
mesh.topology.create_connectivity(2, 0)
mesh.topology.create_connectivity(0, 2)

# Define integration measures
dx = Measure('dx', domain=mesh, subdomain_data=ct, metadata={'quadrature_degree': 2})

# Define function spaces
cell = mesh.ufl_cell().cellname()
CG_V = element("CG", cell, 2, shape=())
DG_0 = element("DG", cell, 0, shape=())
V = functionspace(mesh, CG_V)
V_DG = functionspace(mesh, DG_0)

# Define Dirichlet boundary conditions
def on_dirichlet_bound(x):
    return np.isclose(x[0] ** 2 + x[1] ** 2, (40*1e-3) ** 2)
bcs = []
tdim = mesh.topology.dim - 1
facets = locate_entities_boundary(mesh, dim=tdim, marker=on_dirichlet_bound)
if len(facets) == 0:
    raise ValueError("Dirichlet boundary conditions not defined")
dofs = locate_dofs_topological(V, tdim, facets)
bcs.extend([dirichletbc(default_scalar_type(0.0), dofs, V)])

# Define coils
poles = 8  # Number of poles in the machine
periodicity = 1 # Periodicity of the machine geometry
offset = {"A": 0.0, "B": - 2 / 3 * np.pi, "C": + 2 / 3 * np.pi}
phase_dir = np.tile([-1, 1], round(poles / periodicity) // 2 + 1)[:round(poles / periodicity)]
coils = {f"Phase{p}{i + 1}": {"offset": offset[p] + np.deg2rad(20),
                              "direction": phase_dir[i] if p in "AB" else -phase_dir[i]}
         for i in range(round(poles / periodicity)) for p in "ACB"}
J = Function(V_DG)
J.x.array[:] = 0.0

for i, phase in enumerate(coils.keys()):
    geom = geometry[phase]["geom"]
    cells = ct.find(geom["tag"])

    I_coil = as_ufl(default_scalar_type(coils[phase]["direction"]) * I_ph)  # [A] Current magnitude
    i_coil = I_coil * ufl.sin(position + coils[phase]["offset"])
    J.x.array[cells] = i_coil / assemble_scalar(form(2 / len(geom["surface"]) * dx(geom["tag"])))
# Apply forward scatter to ensure correct parallel distribution
J.x.scatter_forward()

# Define functions to solve
Az = Function(V)
u, v = TrialFunction(V), TestFunction(V)

L = J * v * dx

# Define linear solver
R = nu(Az) * inner(grad(u), grad(v)) * dx
from dolfinx.fem.petsc import LinearProblem
linear_problem = LinearProblem(R, L, u=Az, bcs=bcs,
                               petsc_options={"ksp_type": "gmres",
                                              "ksp_rtol": 1e-6,
                                              "ksp_atol": 1e-10,
                                              "ksp_max_it": 1000})
linear_problem.solve()
new_nu = nu(Az) # Update iron reluctivity with solution of the linear problem

# Define nonlinear solver
R = new_nu * inner(curl(Az), curl(v)) * dx
F = R - L
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
problem = NonlinearProblem(F, u=Az, bcs=bcs)
solver = NewtonSolver(mesh.comm, problem)

opts = PETSc.Options() # type: ignore
ksp = solver.krylov_solver
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "gmres"
opts[f"{option_prefix}ksp_rtol"] = 1.0e-6
opts[f"{option_prefix}ksp_atol"] = 1.0e-8
opts[f"{option_prefix}ksp_min_it"] = 5
opts[f"{option_prefix}ksp_max_it"] = 500
opts[f"{option_prefix}ksp_reuse_preconditioner"] = True
opts[f"{option_prefix}ksp_monitor"] = ""  # Enable the default monitor to stdout
opts[f"{option_prefix}ksp_error_if_not_converged"] = True
opts[f"{option_prefix}ksp_converged_rate"] = ""

# Preconditioner settings
opts[f"{option_prefix}pc_type"] = "gamg"
opts[f"{option_prefix}pc_gamg_threshold"] = 0.5
ksp.setFromOptions()

excitation_fractions = np.linspace(0, 1, 5)[1:]
try:
    for f in excitation_fractions:
        L = f * J * v * dx
        F = R - L
        n, converged = solver.solve(Az)
        print(f"Number of iterations: {n:d}")
        nu_new = nu(Az)
except RuntimeError:
    print("Failed to converge!")
    quit()


## ------ Plot output ------
# Define cmap of Ansys Maxwell
colors = ["ff0000", "ff4900", "ff9200", "ffdb00", "dbff00", "92ff00", "49ff00", "00ff00",
          "00ff49", "00ff92", "00ffdb", "00dbff", "0092ff", "0049ff", "0000ff",]
def hex_to_rgba(hex_color):
    hex_color = hex_color.lstrip('#')
    if len(hex_color) == 6:
        r = int(hex_color[0:2], 16) / 255.0
        g = int(hex_color[2:4], 16) / 255.0
        b = int(hex_color[4:6], 16) / 255.0
        return r, g, b, 1.0
    return 0, 0, 0, 1.0
colors = [hex_to_rgba(c) for c in colors[::-1]]
amc = matplotlib.colors.ListedColormap(colors, "Ansys maxwell colormap")

# Plot current density
plotter = pyvista.Plotter()
grid = pyvista.UnstructuredGrid(*vtk_mesh(mesh, mesh.topology.dim))
grid.cell_data["Jz"] = J.x.array
grid.set_active_scalars("Jz")

Jz_max = np.max([abs(J.x.array.min()), abs(J.x.array.max())])
plotter.add_mesh(grid, show_edges=False, cmap=amc, clim=[-Jz_max, Jz_max])
plotter.view_xy()
plotter.show()
plotter.screenshot("Jz.png")

# Plot vector potential
plotter = pyvista.Plotter()
Az_grid = pyvista.UnstructuredGrid(*vtk_mesh(V))
Az_grid.point_data["A_z"] = Az.x.array
Az_grid.set_active_scalars("A_z")
plotter.add_mesh(Az_grid, show_edges=False, cmap=amc)
plotter.add_mesh(Az_grid.extract_all_edges(), color="grey", opacity=0.25)
plotter.view_xy()
plotter.show()
plotter.screenshot("Az.png")

# Plot flux density
plotter = pyvista.Plotter()
grid = pyvista.UnstructuredGrid(*vtk_mesh(V))
B_func = Function(V)
B_expr = Expression(ufl.sqrt(Az.dx(1) ** 2 + (Az.dx(0)) ** 2) / ufl.sqrt(2), V.element.interpolation_points())
B_func.interpolate(B_expr)

num_dofs = V.dofmap.index_map.size_local
values = np.zeros((num_dofs, 3), dtype=np.float64)
values[:, :V.mesh.geometry.dim] = B_func.x.array.real.reshape(num_dofs, V.dofmap.index_map_bs)

grid = pyvista.UnstructuredGrid(*vtk_mesh(V))
grid.point_data["B"] = values
grid.set_active_vectors("B")
plotter.add_mesh(grid, show_edges=False, cmap=amc, clim=[0, 1.5])
plotter.add_mesh(Az_grid.extract_all_edges(), color="grey", opacity=0.25)
plotter.view_xy()
plotter.show()
plotter.screenshot("B.png")

# Plot permeability
num_cells = V.mesh.topology.index_map(V.mesh.topology.dim).size_local
nu_hd = Function(V_DG)
nu_hd.interpolate(nu(Az))
cell_values = 1 / (nu_hd.x.array.real.reshape((num_cells, V.dofmap.index_map_bs)) * 4*np.pi*1e-7)

grid = pyvista.UnstructuredGrid(*vtk_mesh(V.mesh, V.mesh.topology.dim))
grid.cell_data["μ_r"] = cell_values.flatten()
grid.set_active_scalars("μ_r")

plotter = pyvista.Plotter()
plotter.add_mesh(grid, show_edges=False, cmap=amc, clim=[0, 7000])
plotter.add_mesh(Az_grid.extract_all_edges(), color="grey", opacity=0.25)
plotter.view_xy()
plotter.show()
plotter.screenshot("mu.png")
5

Thanks for the reference! I’m gonna have a look!

6

In dolfiny, they implemented a physical unit system. Maybe this is of use to you?

See these slides:

And here is a demo: dolfiny/demo/units/navier_stokes.py at main · fenics-dolfiny/dolfiny · GitHub

7

In an effort to get more feedback from the nonlinear solver I implemented a custom newton solver based on Custom Newton solvers — FEniCSx tutorial. This allows me to update the constitutive relation (nu) very strictly albeit not very elegantly.

import pickle

import matplotlib
import numpy as np
from numpy import pi
import pyvista
import ufl
from basix.ufl import element
from dolfinx import default_scalar_type
from dolfinx.fem import functionspace, locate_dofs_topological, dirichletbc, Function, Expression, form, assemble_scalar
import dolfinx.fem.petsc
from dolfinx.io.gmshio import read_from_msh
from dolfinx.mesh import locate_entities_boundary
from dolfinx.plot import vtk_mesh
from mpi4py.MPI import COMM_WORLD
from petsc4py import PETSc
from scipy.interpolate import PchipInterpolator
from ufl import Measure, TrialFunction, TestFunction, as_ufl, inner, curl, grad, as_vector

I_ph = 500*9
position = 0.0

bh_curve =  [(0.0000, 0.0000), (0.2447, 38.9522), (0.5008, 60.7621), (0.7509, 104.1758), (1.0000, 187.0970), (1.2522, 426.7090),
            (1.3949, 1157.2412), (1.4994, 2804.8711), (1.5934, 5220.9581), (1.7116, 10592.0), (1.7507, 13305.0),
            (1.7917, 16906.0), (1.8303, 21753.0), (1.8486, 25194.0), (1.8664, 29701.0), (1.8883, 37281.0),
            (1.9057, 45251.0), (1.9230, 55000.0), (1.9404, 66226.0), (1.9578, 78595.0), (1.9751, 91776.0),
            (2.0055, 115790.0), (2.1117, 200310.0), (2.4947, 505020.0), (3.5124, 1314900.0)]

def nu(u: Function) -> Function:
    from scipy.interpolate import CubicSpline
    mu0 = 4 * np.pi * 1e-7  # Permeability of free space
    B_vals, H_vals = zip(*bh_curve)
    B_vals = np.array([b for b in B_vals])
    H_vals = np.array([h for h in H_vals])
    y = np.diff(np.insert(H_vals, 0, -H_vals[1])) / np.diff(np.insert(B_vals, 0, -B_vals[1]))
    y = np.clip(y, a_min=1/7550 * 1/mu0, a_max=None)

    # cs = CubicSpline(B_vals, y, bc_type='clamped')
    cs = PchipInterpolator(B_vals, y)
    B_breaks, B_coefs = cs.x, cs.c

    dm = V_DG.dofmap
    c2d = dm.cell_dofs

    nu_f = Function(V_DG)
    nu_f.x.array[:] = 1 / mu0
    nu_func = Function(V_DG)
    nu_func.x.array[:] = 1 / mu0

    B_norm = ufl.sqrt(u.dx(1) ** 2 + (-u.dx(0)) ** 2)  / ufl.sqrt(2)

    nu_interp = 0
    for i in range(len(B_breaks) - 1):
        # mask = (B_norm >= B_breaks[i]) & (B_norm < B_breaks[i + 1])
        delta_x = B_norm - B_breaks[i]
        if i == len(B_breaks) - 2:
            # Completely saturated
            poly = 1 / mu0
            mask = B_norm >= B_breaks[i]
        else:
            poly = (B_coefs[0, i] * delta_x ** 3 +
                    B_coefs[1, i] * delta_x ** 2 +
                    B_coefs[2, i] * delta_x +
                    B_coefs[3, i])
            mask = ufl.And(B_norm >= B_breaks[i], B_norm < B_breaks[i + 1])
        nu_interp += ufl.conditional(mask, poly, 0)
    nu_expr = Expression(nu_interp, V_DG.element.interpolation_points())
    nu_func.interpolate(nu_expr)

    for key in geometry.keys() - {"Center"}:
        out = geometry[key]
        tag = out["geom"]["tag"]

        if key in ["Stator", "Rotor"]:
            for c in ct.find(tag):
                nu_f.x.array[c2d(c)] = nu_func.x.array[c2d(c)]
        else:
            mur = out["material"]["mu_r"]
            nu_val = 1 / (mu0 * mur)
            for c in ct.find(tag):
                nu_f.x.array[c2d(c)] = nu_val
    nu_f.x.scatter_forward()
    return nu_f

# Load the geometry pickle
with open("geometry.pkl", 'rb') as f:
    geometry = pickle.load(f)
# mesh, ct, ft = read_from_msh("full_geometry.msh", comm=COMM_WORLD, rank=0, gdim=2)
mesh, ct, ft = read_from_msh("simple_ipm_machine.msh", comm=COMM_WORLD, rank=0, gdim=2)
mesh.topology.create_connectivity(2, 1)
mesh.topology.create_connectivity(1, 2)
mesh.topology.create_connectivity(2, 0)
mesh.topology.create_connectivity(0, 2)

# Define integration measures
dx = Measure('dx', domain=mesh, subdomain_data=ct, metadata={'quadrature_degree': 2})

# Define function spaces
cell = mesh.ufl_cell().cellname()
CG_V = element("CG", cell, 2, shape=())
DG_0 = element("DG", cell, 0, shape=())
V = functionspace(mesh, CG_V)
V_DG = functionspace(mesh, DG_0)

# Define Dirichlet boundary conditions
def on_dirichlet_bound(x):
    return np.isclose(x[0] ** 2 + x[1] ** 2, (40*1e-3) ** 2)
bcs = []
tdim = mesh.topology.dim - 1
facets = locate_entities_boundary(mesh, dim=tdim, marker=on_dirichlet_bound)
if len(facets) == 0:
    raise ValueError("Dirichlet boundary conditions not defined")
dofs = locate_dofs_topological(V, tdim, facets)
bcs.extend([dirichletbc(default_scalar_type(0.0), dofs, V)])

# Define coils
poles = 8  # Number of poles in the machine
periodicity = 1 # Periodicity of the machine geometry
offset = {"A": 0.0, "B": - 2 / 3 * np.pi, "C": + 2 / 3 * np.pi}
phase_dir = np.tile([-1, 1], round(poles / periodicity) // 2 + 1)[:round(poles / periodicity)]
coils = {f"Phase{p}{i + 1}": {"offset": offset[p] + np.deg2rad(20),
                              "direction": phase_dir[i] if p in "AB" else -phase_dir[i]}
         for i in range(round(poles / periodicity)) for p in "ACB"}
Jz = Function(V_DG)
Jz.x.array[:] = 0.0

for i, phase in enumerate(coils.keys()):
    geom = geometry[phase]["geom"]
    cells = ct.find(geom["tag"])

    I_coil = as_ufl(default_scalar_type(coils[phase]["direction"]) * I_ph)  # [A] Current magnitude
    i_coil = I_coil * ufl.sin(position + coils[phase]["offset"])
    Jz.x.array[cells] = i_coil / assemble_scalar(form(2 / len(geom["surface"]) * dx(geom["tag"])))
# Apply forward scatter to ensure correct parallel distribution
Jz.x.scatter_forward()

# Define magnets
# Material properties
mu0 = 4 * pi * 1e-7
H_c: float = -920e3  # [A/m] Coercive force
B_r: float = 1.2  # [T] Remanent flux density
B_r_parallel: float = 1.16  # [T] Remanent flux density parallel to the magnetization
B_r_perpendicular: float = 1.22  # [T] Remanent

# Magnet order & Angle
list1 = ["Magnet_1", "Magnet_3", "Magnet_5", "Magnet_7"]
list2 = ["Magnet_2", "Magnet_4", "Magnet_6", "Magnet_8"]
magnet_order = [x for pair in zip(list1, list2) for x in pair]
list1 = [pi + 2 * pi / 8 * i for i in range(0, 8, 2)]
list2 = [2 * pi / 8 * i for i in range(1, 8, 2)]
magnet_angle = [x for pair in zip(list1, list2) for x in pair]

Mx = Function(V_DG)
My = Function(V_DG)
Mx.x.array[:] = 0.0
My.x.array[:] = 0.0

# Loop over all magnets and set the magnetization to the correct vector
for i, pm in enumerate(magnet_order):
    if pm not in geometry.keys():
        raise ValueError("Magnet key not defined")

    cells = ct.find(geometry[pm]["geom"]["tag"])

    vector = as_vector([ufl.cos(magnet_angle[i]), ufl.sin(magnet_angle[i])])
    M_pm_parallel = (B_r_parallel / mu0) * vector[0]
    M_pm_perpendicular = (B_r_perpendicular / mu0) * vector[1]
    Mx.x.array[cells] = M_pm_parallel
    My.x.array[cells] = M_pm_perpendicular
Mx.x.scatter_forward()
My.x.scatter_forward()
M_pm = as_vector([Mx, My])

# Define functions to solve
Az = Function(V)
u, v = TrialFunction(V), TestFunction(V)

n_factor = 1 # Normalization factor
L = n_factor * (Jz * v * dx + mu0 * nu(Az) * inner(M_pm, curl(v)) * dx)

# Define linear solver
R = 1/n_factor * nu(Az) * inner(grad(u), grad(v)) * dx
from dolfinx.fem.petsc import LinearProblem
linear_problem = LinearProblem(R, L, u=Az, bcs=bcs,
                               petsc_options={"ksp_type": "gmres",
                                              "ksp_rtol": 1e-6,
                                              "ksp_atol": 1e-10,
                                              "ksp_max_it": 1000})
linear_problem.solve()

nu_func = Function(V_DG)
nu_func.x.array[:] = nu(Az).x.array[:]
nu_func.x.scatter_forward()

# Define nonlinear solver
R = n_factor * nu_func * inner(curl(Az), curl(v)) * dx
L = n_factor * (Jz * v * dx + mu0 * nu_func * inner(M_pm, curl(v)) * dx)
F = R - L
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
problem = NonlinearProblem(F, u=Az, bcs=bcs)

solver_type = ["default", "custom"][0]

if solver_type == "default":
    solver = NewtonSolver(mesh.comm, problem)
    solver.atol = 1e3
    solver.rtol = 1e3

    opts = PETSc.Options() # type: ignorex
    ksp = solver.krylov_solver
    option_prefix = ksp.getOptionsPrefix()
    opts[f"{option_prefix}ksp_type"] = "gmres"
    opts[f"{option_prefix}ksp_rtol"] = 1.0e3
    opts[f"{option_prefix}ksp_atol"] = 1.0e3
    opts[f"{option_prefix}ksp_min_it"] = 5
    opts[f"{option_prefix}ksp_max_it"] = 500
    opts[f"{option_prefix}ksp_reuse_preconditioner"] = True
    # opts[f"{option_prefix}ksp_monitor"] = ""  # Enable the default monitor to stdout
    opts[f"{option_prefix}ksp_error_if_not_converged"] = True
    # opts[f"{option_prefix}ksp_converged_rate"] = ""

    # Preconditioner settings
    opts[f"{option_prefix}pc_type"] = "gamg"
    opts[f"{option_prefix}pc_gamg_threshold"] = 0.5
    ksp.setFromOptions()

    excitation_fractions = np.linspace(0, 1, 5)[1:]
    try:
        for f in excitation_fractions:
            nu_func.x.array[:] = nu(Az).x.array[:]
            nu_func.x.scatter_forward()
            L = n_factor * (Jz * v * dx + mu0 * nu_func * inner(M_pm, curl(v)) * dx)
            F = R - L
            n, converged = solver.solve(Az)
            print(f"Number of iterations: {n:d}")
            nu_new = nu(Az)
    except RuntimeError:
        print("Failed to converge!")
        quit()
else:
    V = Az.function_space
    residual = dolfinx.fem.form(F)
    J = ufl.derivative(F, Az)
    jacobian = dolfinx.fem.form(J)

    A = dolfinx.fem.petsc.create_matrix(jacobian)
    L = dolfinx.fem.petsc.create_vector(residual)

    solver = PETSc.KSP().create(mesh.comm)
    solver.setOperators(A)
    du = dolfinx.fem.Function(V)

    i = 0
    du_norm = []

    max_iterations = 25
    coords = V.tabulate_dof_coordinates()[:, 0]
    sort_order = np.argsort(coords)
    solutions = np.zeros((max_iterations + 1, len(coords)))
    solutions[0] = Az.x.array[sort_order]

    while i < max_iterations:
        # Assemble Jacobian and residual
        with L.localForm() as loc_L:
            loc_L.set(0)
        A.zeroEntries()
        dolfinx.fem.petsc.assemble_matrix(A, jacobian, bcs=bcs)
        A.assemble()
        dolfinx.fem.petsc.assemble_vector(L, residual)
        L.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)

        # Scale residual by -1
        L.scale(-1)

        # Compute b - J(u_D-u(i-1))
        dolfinx.fem.petsc.apply_lifting(L, [jacobian], [bcs], x0=[Az.x.petsc_vec], alpha=1)
        # Set du|_bc = u_{i-1}-u_D
        dolfinx.fem.petsc.set_bc(L, bcs, Az.x.petsc_vec, 1.0)
        L.ghostUpdate(addv=PETSc.InsertMode.INSERT_VALUES, mode=PETSc.ScatterMode.FORWARD)

        # Solve linear problem
        solver.solve(L, du.x.petsc_vec)
        du.x.scatter_forward()

        # Update u_{i+1} = u_i + delta u_i
        Az.x.array[:] += 0.8 ** i * du.x.array
        i += 1

        # Update reluctivity
        nu_func.x.array[:] = nu(Az).x.array[:]
        nu_func.x.scatter_forward()

        # Compute norm of update
        correction_norm = du.x.petsc_vec.norm(0)
        du_norm.append(correction_norm)
        print(f"Iteration {i}: Correction norm {correction_norm}")
        if correction_norm < 1e-10:
            break
        solutions[i, :] = Az.x.array[sort_order]

    dolfinx.fem.petsc.assemble_vector(L, residual)
    print(f"Final residual {L.norm(0)}")
    A.destroy()
    L.destroy()
    solver.destroy()


Az.x.array[:] = Az.x.array[:] / n_factor
## ------ Plot output ------
# Define cmap of Ansys Maxwell
colors = ["ff0000", "ff4900", "ff9200", "ffdb00", "dbff00", "92ff00", "49ff00", "00ff00",
          "00ff49", "00ff92", "00ffdb", "00dbff", "0092ff", "0049ff", "0000ff",]
def hex_to_rgba(hex_color):
    hex_color = hex_color.lstrip('#')
    if len(hex_color) == 6:
        r = int(hex_color[0:2], 16) / 255.0
        g = int(hex_color[2:4], 16) / 255.0
        b = int(hex_color[4:6], 16) / 255.0
        return r, g, b, 1.0
    return 0, 0, 0, 1.0
colors = [hex_to_rgba(c) for c in colors[::-1]]
amc = matplotlib.colors.ListedColormap(colors, "Ansys maxwell colormap")

# Plot current density
plotter = pyvista.Plotter()
grid = pyvista.UnstructuredGrid(*vtk_mesh(mesh, mesh.topology.dim))
grid.cell_data["Jz"] = Jz.x.array
grid.set_active_scalars("Jz")

Jz_max = np.max([abs(Jz.x.array.min()), abs(Jz.x.array.max())])
plotter.add_mesh(grid, show_edges=False, cmap=amc, clim=[-Jz_max, Jz_max])
plotter.view_xy()
plotter.show()
plotter.screenshot("Jz.png")

# Plot vector potential
plotter = pyvista.Plotter()
Az_grid = pyvista.UnstructuredGrid(*vtk_mesh(V))
Az_grid.point_data["A_z"] = Az.x.array
Az_grid.set_active_scalars("A_z")
plotter.add_mesh(Az_grid, show_edges=False, cmap=amc)
plotter.add_mesh(Az_grid.extract_all_edges(), color="grey", opacity=0.25)
plotter.view_xy()
plotter.show()
plotter.screenshot("Az.png")

# Plot flux density
plotter = pyvista.Plotter()
grid = pyvista.UnstructuredGrid(*vtk_mesh(V))
B_func = Function(V)
B_expr = Expression(ufl.sqrt(Az.dx(1) ** 2 + (Az.dx(0)) ** 2) / ufl.sqrt(2), V.element.interpolation_points())
B_func.interpolate(B_expr)

num_dofs = V.dofmap.index_map.size_local
values = np.zeros((num_dofs, 3), dtype=np.float64)
values[:, :V.mesh.geometry.dim] = B_func.x.array.real.reshape(num_dofs, V.dofmap.index_map_bs)

grid = pyvista.UnstructuredGrid(*vtk_mesh(V))
grid.point_data["B"] = values
grid.set_active_vectors("B")
plotter.add_mesh(grid, show_edges=False, cmap=amc, clim=[0, 1.5])
plotter.add_mesh(Az_grid.extract_all_edges(), color="grey", opacity=0.25)
plotter.view_xy()
plotter.show()
plotter.screenshot("B.png")

# Plot permeability
num_cells = V.mesh.topology.index_map(V.mesh.topology.dim).size_local
nu_hd = Function(V_DG)
nu_hd.interpolate(nu(Az))
cell_values = 1 / (nu_hd.x.array.real.reshape((num_cells, V.dofmap.index_map_bs)) * 4*np.pi*1e-7)

grid = pyvista.UnstructuredGrid(*vtk_mesh(V.mesh, V.mesh.topology.dim))
grid.cell_data["μ_r"] = cell_values.flatten()
grid.set_active_scalars("μ_r")

plotter = pyvista.Plotter()
plotter.add_mesh(grid, show_edges=False, cmap=amc, clim=[0, 7000])
plotter.add_mesh(Az_grid.extract_all_edges(), color="grey", opacity=0.25)
plotter.view_xy()
plotter.show()
plotter.screenshot("mu.png")

The result of this is that the custom solver quite rapidly converges to some close approximation of the solution, with the saturation effect clearly visible in early results. Compared to the default NewtonSolver which doesn’t converge at all.

Clearly the custom solver would need a lot of refinement to efficiently solve, which is why I would prefer using a SNES solver. The custom solver would need to update this nu_func function during solving and sent some information on the residual to the console. Can this be achieved by building a CustomNewtonSolver class? Or perhaps by interfacing with PETSc directly?
Any insight or sources would be appreciated.

2 Likes
8

You could overload the function passed to: snes.setFunction and snes.setJacobian to call your update, ref:

i.e. something like

import functools
def custom_assemble_residual(u, residual, jacobian, bcs, _snes, x,b):
    nu_func.x.array[:] = nu(Az).x.array[:]
    nu_func.x.scatter_forward()
    return dolfinx.fem.petsc.assemble_residual(u, residual, jacobian, bcs, _snes, x, b)

problem = dolfinx.fem.petsc.NonlinearProblem(...)
problem.solver.setFunction(functools.partial(custom_assemble_residual, problem.u, problem.F, problem.J, bcs), problem.b)
1 Like
9

Thanks dokken,

That’s exactly what I am trying to achieve for my problem. I would like to add a stabilization term, (which I cannot express as a ufl form.), to the system matrix (Jacobian) of my nonlinear problem.

But it seems like this recipe is depending on the implementation details of this version of dolfinx. Hope that it provides public hooks for user customization.

:folded_hands:

10

As the functions assemble_jacobian and assemble_residual (dolfinx/python/dolfinx/fem/petsc.py at main · FEniCS/dolfinx · GitHub and
dolfinx/python/dolfinx/fem/petsc.py at main · FEniCS/dolfinx · GitHub)
simply fills in data in J and b respectively, the overloaded function can modify these in a non-ufl manner (as you have access to the vector/matrix etc.)

This interface was introduced in v0.10.0 (i.e. the streamlined interface), but it can be achieved in older versions of FEniCSx as well by implementing the assemble_jacobian and assemble_residual in that version of the code.

11

@dokken
I think I also need to modify self._A = create_matrix(self.J, kind=kind) since the modified Jacobian will have a denser sparsity pattern.

12

Thanks Dokken!

Once I’ve managed to switch to dolfinx 0.10.0 I’ll give this a try.

I think the new interface will make this significantly easier.

13

@JvLuenen
I didn’t run your code. But it seems like you use Lagrange elements for ufl.curl. Do you have any specific reason to choose the Lagrange elements instead of N1curl elements. I guess your geometry is regular enough, it would not cause any issues. Do your geometry has sharp corners larger then 180 degrees?

14

@young What would be the reason to use N1curl elements? I don’t have any sharp corners > 180 degrees so maybe it would be possible. Currently changing the element type to “N1curl” But this just raises the error:

’’’

File “/dolfinx-env/lib/python3.12/site-packages/basix/ufl.py”, line 1915, in blocked_element
raise ValueError(“Cannot create a blocked element containing a non-scalar element.”)
ValueError: Cannot create a blocked element containing a non-scalar element.

‘‘‘

15

In 2D geometry for in-plane magnetic field, the problem becomes scalar in terms of the vector potential. You can verify it on paper by evaluating B=\nabla\times A with derivatives out of the plane to be zero. Moreover, A_z is continuous unlike B. That’s why using Lagrange elements are adequate for approximating A_z

1 Like
16

@dokken I have implemented your proposed solution to update the residual and did so for the jacobian as well. The snes newtonls solver with preonly ksp solver does not manage to solve the problem. I think this is caused by the large drop-off in the constitutive relationship.

SNES Object: (magnetostatics_nonlinear_) 1 MPI process
  type: newtonls
  SNES has not been set up so information may be incomplete
  maximum iterations=25, maximum function evaluations=10000
  tolerances: relative=1e-08, absolute=1e-08, solution=1e-08
  total number of linear solver iterations=0
  total number of function evaluations=0
  norm schedule ALWAYS
  SNESLineSearch Object: (magnetostatics_nonlinear_) 1 MPI process
    type: basic
    maxlambda=1.000000e+00, minlambda=1.000000e-12
    tolerances: relative=1.000000e-08, absolute=1.000000e-15, lambda=1.000000e-08
    maximum iterations=40
  KSP Object: (magnetostatics_nonlinear_) 1 MPI process
    type: preonly
    maximum iterations=10000, initial guess is zero
    tolerances: relative=1e-05, absolute=1e-50, divergence=10000.
    left preconditioning
    using DEFAULT norm type for convergence test
  PC Object: (magnetostatics_nonlinear_) 1 MPI process
    type: lu
    PC has not been set up so information may be incomplete
      out-of-place factorization
      tolerance for zero pivot 2.22045e-14
    linear system matrix = precond matrix:
    Matrix has not been assembled yet
  0 SNES Function norm 1.209358416584e+04
  1 SNES Function norm 2.418716833167e+03
  2 SNES Function norm 2.689237647745e+05
  3 SNES Function norm 5.415757832413e+04
  4 SNES Function norm 4.177713053994e+04
  5 SNES Function norm 5.423290235571e+07
  6 SNES Function norm 1.084658234680e+07
  7 SNES Function norm 2.170081025475e+06
  8 SNES Function norm 2.346019747068e+07
  9 SNES Function norm 8.954619164968e+07
 10 SNES Function norm 9.334662983160e+07
 11 SNES Function norm 1.609544927951e+08
  Nonlinear magnetostatics_nonlinear_ solve did not converge due to DIVERGED_DTOL iterations 11

Forcing the newton steps to be very small, by setting the damping to 0.05 does push the solver towards the right solution. Do you have any ideas on some solver settings what would be able to cope with this highly nonlinear relationship?