Hi everyone,
I’ve been struggling with implementing nonlinear materials into maxwells equations using the mpc nonlinear problem. To simplify the problem I have removed periodicity, allowing me to use the doflinx.fem.petsc.NonlinearProblem.
As the results were still not as expected I removed the nonlinear material property and assumed it to be linear for now. I also removed any MMF source from the problem (so no coils or magnets), and instead set a dirichlet boundary condition to a 0.01 for the magnetic vector potential that I am solving for. Comparing the linear solver, and nonlinear solver the nonlinear solver arrives at a vastly different result however.
The linear solution is correct, given that the disc has a permeability of 3000 with slots of permeability=1, reducing the magnetic vector potential inside the material. The nonlinear solution is practically a homogeneous field, which is incorrect.
Below is a MWE
To create the mesh:
from dataclasses import dataclass
from math import pi, radians, cos, sin, tan, ceil
import gmsh
from sympy import symbols, Eq, solve
mm = 1e-3 # [m] millimeter to meter conversion factor
@dataclass
class MachineDims:
dimension: int = 2 # [int] Dimension of the geometry, 2D or 3D
periodicity: int = 1
geom_start_angle: float = 0.0 # [rad] Start angle of the geometry
geom_end_angle: float = 1/periodicity * 2 * pi # [rad] end angle of the geometry in radians
l_ag: float = 1.5*mm # [m] Radial length of the air gap
slots = 48 # number of slots
@dataclass
class WedgeDims:
r_inner: float
r_outer: float
@dataclass
class StatorDims(WedgeDims):
r_inner: float = 161.90/2*mm
r_outer: float = 269.24/2*mm
@dataclass
class RotorDims(WedgeDims):
r_inner: float = 110.64/2*mm
r_outer: float = 160.40/2*mm
@dataclass
class SlotDims:
a_slot_gap = radians(1.366) # [rad]
h_tooth_tip = 1.03*mm # [m] height of the tooth tip
w_tooth_tip = 1.535*mm # [m] width of the tooth tip
h_slot = 29.50*mm # [m] height of the slot
a_slot_widening = radians(1.25) # [rad] angle of the slot widening
h_insulation = 3.50*mm # [m] height of the insulation
@dataclass
class DuctDims:
a_duct = 65.0 / 180 * pi # [rad]
r_duct_maximum = 78.72*mm # [m]
r_duct_minimum = 62.6*mm # [m]
h_duct = 4.7*mm # [m]
o_magnet = 2.25*mm # [m] offset of the magnet from the mirror plane
h_magnet = 6.48*mm # [m] height of the magnet
w_magnet = 16.0*mm # [m] width of the magnet
w_duct_peak = 3.0*mm # [m] width of the duct peak
def mirror_entities(entities, angle, point=(0,0,0)):
"""
Mirror a list of entities across a plane defined by a normal vector and a point.
:param entities: List of entities to be mirrored.
:param angle: Angle to the plane to the x-axis.
:param point: A point on the plane across which to mirror.
:return: List of mirrored entities.
"""
a = sin(angle)
b = -cos(angle)
c = 0
d = -(a * point[0] + b * point[1] + c * point[2])
e = factory.copy(entities)
factory.mirror(e, a, b, c, d)
return e
def repeat_shape_on_angle(out: list[tuple], occurrence: int = None, periodicity: int = None, angle_offset: float = None,
angle: float = None) -> list:
"""
Repeat a shape defined in out[0] around the z-axis at an angle defined by the number of slots and periodicity.
:param out: List containing the shape to be repeated.
:param occurrence: Total number of slots in the machine.
:param periodicity: Number of slots in one period.
:param angle_offset: Offset angle to apply to each repetition.
:param angle: Specific angle to rotate the shape. If None, it will be calculated based on slots and periodicity.
:return: List of entities after rotation.
"""
output = []
for k, (order, shape) in enumerate(out):
if order != 2:
raise ValueError(f"Expected a surface (2), but got {order} at index {k}.")
if occurrence is not None and periodicity is not None:
out_new = [shape]
for i in range(1, ceil(occurrence / periodicity)):
if angle_offset is None:
angle_offset = 0
angle = i * (2 * pi / occurrence) + angle_offset
rotated_surface = factory.copy([(order, shape)])
factory.rotate(rotated_surface, 0, 0, 0, 0, 0, 1, angle)
out_new.append(rotated_surface[0][1])
elif angle is not None:
rotated_surface = factory.copy([(order, shape)])
factory.rotate(rotated_surface, 0, 0, 0, 0, 0, 1, angle)
out_new = [rotated_surface[0][1]]
else:
raise ValueError("Either angle or both slots and periodicity must be provided.")
output.extend(out_new)
return output
def cut_from_surface(out: list, tool_surfaces: list) -> list:
"""
Cut the tool surfaces from the output surface.
:param out: List containing the surface to be cut.
:param tool_surfaces: List of surfaces to be used as tools for cutting.
:return: List of entities after cutting.
"""
surfaces_to_cut = [(2, surface) for surface in tool_surfaces]
cut_result, _ = factory.cut([(2, out[0])], surfaces_to_cut, removeObject=False, removeTool=False)
return [tag for dim, tag in cut_result if dim == 2] # Keep only the surface tags
def square_in_polar_coordinates(r1, r2, a1, a2) -> tuple:
"""
Create a square in polar coordinates defined by two radii and two angles.
:param r1: Inner radius of the square.
:param r2: Outer radius of the square.
:param a1: Start angle of the square.
:param a2: End angle of the square.
:return: List of points defining the square in polar coordinates.
"""
x = [r1 * cos(a1), r1 * cos(a2), r2 * cos(a2), r2 * cos(a1)]
y = [r1 * sin(a1), r1 * sin(a2), r2 * sin(a2), r2 * sin(a1)]
z = [0, 0, 0, 0] # All points are in the same plane (z=0)
return x, y, z
def draw_square_from_points(out: dict) -> dict:
"""
Draw a square from the points defined in the out dictionary.
:param out: dictionary containing the points of the square.
:return: dictionary containing the points, curves, curve loop, and surface of the square.
"""
out["curve"].extend([
factory.addLine(out["point"][0], out["point"][1]), # Bottom line
factory.addLine(out["point"][1], out["point"][2]), # Right line
factory.addLine(out["point"][2], out["point"][3]), # Top line
factory.addLine(out["point"][3], out["point"][0]) # Left Line
])
# Create the curve loop and surface for the magnet
out["curve_loop"] = factory.addCurveLoop(out["curve"])
out["surface"] = [factory.addPlaneSurface([out["curve_loop"]])]
return out
def solve_algebraic_system(param, eq1, eq2) -> tuple[float, float]:
"""
Solve a system of algebraic equations.
:param param: The parameters to solve for, typically a tuple of symbols.
:param eq1: First equation in the system.
:param eq2: Second equation in the system.
:return: A tuple containing the solutions for x and y.
"""
solutions = solve((eq1, eq2), param)
if solutions[0][0] > 0 and solutions[0][1] > 0:
return solutions[0]
else:
return solutions[1]
gmsh.initialize()
gmsh.model.add("GeometryBuilder")
machine_dims = MachineDims()
stator_dims = StatorDims()
rotor_dims = RotorDims()
slot_dims = SlotDims()
duct_dims = DuctDims()
## Build the geometry
factory = gmsh.model.occ
out = {}
# Create the center point
center_tag = factory.addPoint(0, 0, 0, 1e-2)
out['center'] = center_tag
coil_out = {"point": [], "curve": [], "curve_loop": None, "surface": []}
slots = machine_dims.slots # Number of slots in the machine
a_slot_gap = slot_dims.a_slot_gap # Angle of the slot gap
h_tooth_tip = slot_dims.h_tooth_tip # Height of the tooth tip
h_slot = slot_dims.h_slot # Height of the slot
h_insulation = slot_dims.h_insulation # Height of the insulation
r1 = stator_dims.r_inner + h_tooth_tip + h_insulation
r2 = stator_dims.r_inner + h_tooth_tip + h_slot
a1 = machine_dims.geom_start_angle + (2*pi / slots) / 2 - a_slot_gap
a2 = machine_dims.geom_start_angle + (2*pi / slots) / 2 + a_slot_gap
# Calculate the coordinates of the coil points
x, y, z = square_in_polar_coordinates(r1, r2, a1, a2)
for i in range(len(x)):
point_tag = factory.addPoint(x[i], y[i], z[i], 1e-3)
coil_out["point"].append(point_tag)
coil_out = draw_square_from_points(coil_out)
coils = {"PhaseA1": coil_out,
"PhaseA2": {"surface": []},
"PhaseC1": {"surface": []},
"PhaseC2": {"surface": []},
"PhaseB1": {"surface": []},
"PhaseB2": {"surface": []}}
# Repeat the coil surface for all slots
for i, coil in enumerate(coils.keys()):
if coil == "PhaseA1":
coils[coil]["name"] = coil # type: ignore
coils[coil]["surface"] = repeat_shape_on_angle([(2, coil_out["surface"][0])],
occurrence=round(slots / (len(coils.keys()))),
periodicity=1,
angle_offset=0)
continue
coils[coil]["surface"] = repeat_shape_on_angle([(2, coil_out["surface"][0])], angle=i * (2 * pi / slots))
coils[coil]["surface"] = repeat_shape_on_angle([(2, coils[coil]["surface"][0])],
occurrence=round(slots / (len(coils.keys()))),
periodicity=1)
coils[coil]["name"] = coil # type: ignore
coil_out = coils
slot_out = {"point": [], "curve": [], "curve_loop": None, "surface": []}
slots = machine_dims.slots # Number of slots in the machine
periodicity = machine_dims.periodicity # Periodicity of the machine geometry
a_slot_gap = slot_dims.a_slot_gap # Angle of the slot gap
h_tooth_tip = slot_dims.h_tooth_tip # Height of the tooth tip
w_tooth_tip = slot_dims.w_tooth_tip # Width of the tooth tip
h_slot = slot_dims.h_slot # Height of the slot
a_slot_widening = slot_dims.a_slot_widening # Angle of the slot widening
r1 = stator_dims.r_inner # Inner radius of the stator
r2 = stator_dims.r_inner + h_tooth_tip
r4 = stator_dims.r_inner + h_tooth_tip + h_slot - 1e-4
a1 = machine_dims.geom_start_angle + (2*pi / slots) / 2 - a_slot_gap # Start angle of the slot
a2 = a1 - 1/2*pi
a3 = a2 - a_slot_widening
a4 = machine_dims.geom_start_angle + (2*pi / slots) / 2
# Calculate the coordinates of the slot points
x = [r1*cos(a1), r2*cos(a1), r2*cos(a1) + w_tooth_tip*cos(a2), r2*cos(a1) + w_tooth_tip*cos(a2) - h_slot*sin(a3), r4*cos(a4),]
y = [r1*sin(a1), r2*sin(a1), r2*sin(a1) + w_tooth_tip*sin(a2), r2*sin(a1) + w_tooth_tip*sin(a2) + h_slot*cos(a3), r4*sin(a4),]
z = [0] * len(x)
for i in range(len(x)):
point_tag = factory.addPoint(x[i], y[i], z[i], 1e-3)
slot_out["point"].append(point_tag)
# Repeat the points on the other half of the slot
points = [(0, x) for x in slot_out["point"][0:-1][::-1]]
mirrored_points = mirror_entities(points, angle=a4, point=(0, 0, 0))
slot_out["point"].extend([x[1] for x in mirrored_points])
# Create the curves for the slot
slot_out["curve"].extend([
factory.addLine(slot_out["point"][0], slot_out["point"][1]), # Inner line
factory.addLine(slot_out["point"][1], slot_out["point"][2]), # Tooth tip line
factory.addLine(slot_out["point"][2], slot_out["point"][3]), # Slot wall
factory.addCircleArc(slot_out["point"][3], slot_out["point"][4], slot_out["point"][5]), # Outer arc
factory.addLine(slot_out["point"][5], slot_out["point"][6]), # Slot wall
factory.addLine(slot_out["point"][6], slot_out["point"][7]), # Tooth tip line
factory.addLine(slot_out["point"][7], slot_out["point"][8]), # Inner line
factory.addLine(slot_out["point"][8], slot_out["point"][0])]) # Closing line
# Create the curve loop and surface for the slot
slot_out["curve_loop"] = factory.addCurveLoop(slot_out["curve"])
slot_out["surface"] = [factory.addPlaneSurface([slot_out["curve_loop"]])]
# Cut the coil out of the slot surface
if coil_out is not None:
for key, coil_surface in coil_out.items():
slot_out["surface"].append(cut_from_surface(slot_out["surface"], coil_surface["surface"])[0])
old_surface = slot_out["surface"][0]
factory.remove([(2, old_surface)], recursive=True)
slot_out["surface"] = slot_out["surface"][1:]
# Repeat the slot surface for all slots
slot_out["surface"] = repeat_shape_on_angle([(2, slot_out["surface"][0])], slots, periodicity)
slot_out["name"] = "Slot"
stator_out = {"point": [], "curve": [], "curve_loop": None, "surface": []}
r1 = stator_dims.r_inner # Inner radius of the stator
r2 = stator_dims.r_outer # Outer radius of the stator
a1 = machine_dims.geom_start_angle # Start angle of the geometry
a2 = machine_dims.geom_end_angle # End angle of the geometry
# Create the points for the stator
x, y, z = square_in_polar_coordinates(r1, r2, a1, a2)
for i in range(4):
point_tag = factory.addPoint(x[i], y[i], z[i], 1e-3)
stator_out["point"].append(point_tag)
# Create the curves for the wedge
stator_out["curve"].append(factory.addCircle(0, 0, 0, r1)) # Inner circle
stator_out["curve"].append(factory.addCircle(0, 0, 0, r2)) # Outer circle
curve_loop, surface = [], []
for i in range(2):
curve_loop.append(factory.addCurveLoop([stator_out["curve"][i]]))
surface.append(factory.addPlaneSurface([curve_loop[i]]))
stator_out["surface"] = cut_from_surface([surface[1]], [surface[0]])
surface_to_cut = [slot_out, coil_out]
for surface_tc in surface_to_cut:
if "PhaseA1" in surface_tc.keys():
for key, coil_surface in surface_tc.items():
if "surface" in coil_surface:
stator_out["surface"] = cut_from_surface(stator_out["surface"], coil_surface["surface"])
else:
raise ValueError(f"Coil {key} does not contain a surface to cut from the slot.")
elif surface_tc["name"] in {"Slot"}:
name = "Stator"
stator_out["surface"] = cut_from_surface(stator_out["surface"], surface_tc["surface"])
stator_out["name"] = name
gmsh.model.occ.synchronize()
geometry = {"Center": out,
"Stator": {"geom": stator_out},
"Slot": {"geom": slot_out},
"PhaseA1": {"geom": coils["PhaseA1"]},
"PhaseA2": {"geom": coils["PhaseA2"]},
"PhaseB1": {"geom": coils["PhaseB1"]},
"PhaseB2": {"geom": coils["PhaseB2"]},
"PhaseC1": {"geom": coils["PhaseC1"]},
"PhaseC2": {"geom": coils["PhaseC2"]}}
for key in geometry.keys():
if key in ["Center", "Bounds"]:
continue
out = geometry[key]["geom"]
out["bbox"] = []
for s in out["surface"]:
out["bbox"].append(factory.getBoundingBox(2, s))
gmsh.model.occ.removeAllDuplicates()
gmsh.model.occ.synchronize()
def is_same_bbox(bbox1, bbox2, tol=1e-6):
return all(abs(a - b) < tol for a, b in zip(bbox1, bbox2))
# Find the updated tags
for key in geometry.keys():
if key in ["Center", "Bounds"]:
continue
out = geometry[key]["geom"]
new_s = []
for i, s in enumerate(out["surface"]):
old_bbox = out["bbox"][i] # Assuming all surfaces have the same bounding box
new_tag = None
for (dim, tag) in factory.getEntities(2):
bbox = factory.getBoundingBox(dim, tag)
if is_same_bbox(bbox, old_bbox):
new_tag = tag
break
new_s.append(new_tag)
geometry[key]["geom"]["surface"] = new_s
for key in geometry.keys():
if key in ["Center", "Bounds"]:
continue
out = geometry[key]["geom"]
out["tag"] = gmsh.model.addPhysicalGroup(2, out["surface"])
gmsh.model.setPhysicalName(2, out["tag"], out["name"])
gmsh.model.mesh.generate(2)
gmsh.write("geometry.msh")
gmsh.fltk.run()
gmsh.finalize()
To define and solve the problem:
from math import pi
import numpy as np
import pyvista
import ufl
from petsc4py import PETSc
from basix.ufl import element
from dolfinx import default_scalar_type
from dolfinx.fem import functionspace, Function, locate_dofs_topological, dirichletbc, Expression
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 ufl import Measure, TestFunction, inner, grad, derivative, TrialFunction
geometry = {'Center': {'center': 1},
'PhaseA1': {'geom': {'curve': [1, 2, 3, 4], 'curve_loop': 1, 'name': 'PhaseA1', 'point': [2, 3, 4, 5], 'surface': [1, 2, 3, 4, 5, 6, 7, 8], 'tag': 3}, 'material': {'mu_0': 4*pi * 1e-7, 'mu_r': 1.0}},
'PhaseA2': {'geom': {'name': 'PhaseA2', 'surface': [9, 10, 11, 12, 13, 14, 15, 16], 'tag': 4}, 'material': {'mu_0': 4*pi * 1e-7, 'mu_r': 1.0}},
'PhaseB1': {'geom': {'name': 'PhaseB1', 'surface': [33, 34, 35, 36, 37, 38, 39, 40], 'tag': 5}, 'material': {'mu_0': 4*pi * 1e-7, 'mu_r': 1.0}},
'PhaseB2': {'geom': {'name': 'PhaseB2', 'surface': [41, 42, 43, 44, 45, 46, 47, 48], 'tag': 6}, 'material': {'mu_0': 4*pi * 1e-7, 'mu_r': 1.0}},
'PhaseC1': {'geom': {'name': 'PhaseC1', 'surface': [17, 18, 19, 20, 21, 22, 23, 24], 'tag': 7}, 'material': {'mu_0': 4*pi * 1e-7, 'mu_r': 1.0}},
'PhaseC2': {'geom': {'name': 'PhaseC2', 'surface': [25, 26, 27, 28, 29, 30, 31, 32], 'tag': 8}, 'material': {'mu_0': 4*pi * 1e-7, 'mu_r': 1.0}},
'Slot': {'geom': {'curve': [193, 194, 195, 196, 197, 198, 199, 200], 'curve_loop': 49, 'name': 'Slot', 'point': [194, 195, 196, 197, 198, 199, 200, 201, 202], 'surface': [102, 104, 106, 108, 110, 112, 114, 116, 118, 120, 122, 124, 126, 128, 130, 132, 134, 136, 138, 140, 142, 144, 146, 148, 150, 152, 154, 156, 158, 160, 162, 164, 166, 168, 170, 172, 174, 176, 178, 180, 182, 184, 186, 188, 190, 192, 194, 196], 'tag': 2}, 'material': {'mu_0': 4*pi * 1e-7, 'mu_r': 1.0}},
'Stator': {'geom': {'curve': [765, 766], 'curve_loop': None, 'name': 'Stator', 'point': [767, 768, 769, 770], 'surface': [101], 'tag': 1}, 'material': {'mu_0': 4*pi * 1e-7, 'mu_r': 3000}},
}
mesh, cell_tags, facet_tags = read_from_msh("geometry.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)
dx = Measure("dx", subdomain_data=cell_tags, domain=mesh, metadata={"quadrature_degree": 2})
cell = mesh.ufl_cell().cellname()
# Define function space
CG_V = element("Lagrange", cell, 2, shape=()) # CG: Continuous Galerkin vector potential (Sized for Az only)
DG_0 = element("DG", cell, 0, shape=()) # DG: Discontinuous Galerkin scalar space
V = functionspace(mesh, CG_V)
V_DG = functionspace(mesh, DG_0)
# Boundary conditions
def dirichlet_bounds(x):
return np.isclose(x[0]**2+x[1]**2, (269.24/2/1000)**2)
tdim = mesh.topology.dim - 1
facets = locate_entities_boundary(mesh, dim=tdim, marker=dirichlet_bounds)
dofs = locate_dofs_topological(V, tdim, facets)
bc = dirichletbc(default_scalar_type(0.1), dofs, V)
bcs = [bc]
# Define variational problem
Az = Function(V)
Az.x.array[:] = 0.0
u = TrialFunction(V)
v = TestFunction(V)
nu = Function(V_DG)
dm = V_DG.dofmap
c2d = dm.cell_dofs
for key in geometry.keys() - {"Center"}:
out = geometry[key]
tag = out["geom"]["tag"]
val = 1 / (out["material"]["mu_0"] * out["material"]["mu_r"])
for c in cell_tags.find(tag):
nu.x.array[c2d(c)] = val
nu.x.scatter_forward()
a = inner(nu * grad(u), grad(v)) * dx
Jz = Function(V)
Jz.x.array[:] = 0
L = Jz * v * dx
## Linear variational problem
from dolfinx.fem.petsc import LinearProblem
problem = LinearProblem(a,
L,
u=Az,
bcs=bcs,
petsc_options={"ksp_type": "gmres",
"ksp_rtol": 1e-6,
"ksp_atol": 1e-10,
"ksp_max_it": 1000})
problem.solve()
## Nonlinear variational problem
# R = inner(nu * grad(Az), grad(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)
# solver.convergence_criterion = "residual"
#
# n, converged = solver.solve(Az)
# print(f"Number of iterations: {n:d}")
## ------ Plot results ------
# 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)
plotter.add_mesh(Az_grid.extract_all_edges(), color="grey", opacity=0.25)
plotter.view_xy()
plotter.show()
# 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 + 1e-12), 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)
plotter.add_mesh(Az_grid.extract_all_edges(), color="grey", opacity=0.25)
plotter.view_xy()
plotter.show()
# Plot permeability
num_cells = V.mesh.topology.index_map(V.mesh.topology.dim).size_local
cell_values = 1 / (nu.x.array.real.reshape(num_cells, V.dofmap.index_map_bs) * 4 * pi * 1e-7)
grid = pyvista.UnstructuredGrid(*vtk_mesh(V.mesh, V.mesh.topology.dim))
grid.cell_data["nu"] = cell_values.flatten()
grid.set_active_scalars("nu")
plotter = pyvista.Plotter()
plotter.add_mesh(grid, show_edges=False, cmap="coolwarm")
plotter.add_mesh(Az_grid.extract_all_edges(), color="grey", opacity=0.25)
plotter.view_xy()
plotter.show()
I’m completely lost as to why the nonlinear solver arrives at this (wrong) result. Does anyone have any ideas?
