Nonlinear variational problem takes too long

1

Hello

I’m solving a nonlinear magnetostatic problem. At first, I could solve the problem, but it takes too long comparing the same problem with FEMM (it takes seconds in FEMM and about 6 mintes in Fenics).

Does anyone know why it’s happening?

# mesh parameters
nmesh = 50 # mesh density
tol = 1e-8 # mesh tolerance 

# geometry parameters
# air box
r_air = 80e-3
z_air = 80e-3

air = Rectangle(Point(tol,-z_air), Point(r_air,z_air))

# Coil
thickness = 25.4e-3
width = 2*25.4e-3
r_inner = 0.5*25.4e-3 #coil inner radius

coil = Rectangle(Point(r_inner,-width/2), Point(r_inner+thickness,width/2))

# Core, the air region inside the coil
core = Rectangle(Point(tol,-width/2), Point(r_inner,width/2))

# domain
domain = air + coil + core

# Set subdomain for wire
domain.set_subdomain(1,coil)
domain.set_subdomain(2,core)

# Create mesh
mesh = generate_mesh(domain, nmesh)

# Steel
Bi_steel = np.array([0, 0.221950, 0.245515, 0.344303, 0.375573, 0.454417, 0.627981, 
                     0.670134, 0.861453, 1.075180, 1.241074, 1.423388, 1.656238, 1.686626, 
                     1.813505, 1.964422, 1.979083, 2.012433, 2.021337, 2.033503, 2.050973, 
                     2.052071, 2.191983, 2.197328, 2.240825, 2.309729, 2.327795, 2.435784])
Hi_steel = np.array([tol, 225.366667, 237.316667, 291.793333, 310.450000, 358.730000, 
                     483.890000, 520.136667, 723.673333, 1071.333333, 1570.566667, 
                     2775.500000, 6290.533333, 7049.866667, 12338.666667, 26304.666667, 
                     28581.000000, 36287.000000, 39022.333333, 43292.333333, 50590.000000, 
                     51118.333333, 134313.333333, 138566.666667, 168803.333333, 223476.666667, 
                     237853.333333, 321480.000000])

mi_steel = np.divide(Bi_steel, Hi_steel)
mi_steel[0] = -Bi_steel[1]*((mi_steel[2] - mi_steel[1])/(Bi_steel[2] - Bi_steel[1])) + mi_steel[1]

f2 = interp1d(Bi_steel, mi_steel, kind='quadratic')

# Define function space
V = FunctionSpace(mesh, 'P', 2)

# Define boundary condition
bc = DirichletBC(V, Constant(0), 'on_boundary')

def mi(phi):
    B = project(as_vector((-(1/r)*phi.dx(1), (1/r)*phi.dx(0))))
    norm_B = norm(B)
    return f2(norm_B)

# Define subdomain markers and integration measure
markers = MeshFunction('size_t', mesh, 2, mesh.domains())
dx = Measure('dx', domain=mesh, subdomain_data=markers)

# Define magnetic permeability
class Permeability(UserExpression):
    def __init__(self, markers, **kwargs):
        super().__init__(**kwargs)
        self.markers = markers
        self.mi_0 = mi0            #coil and air
    
    def eval_cell(self, values, x, cell):
        if self.markers[cell.index] == 2:
            values[0] = mi(phi)
        else:
            values[0] = self.mi_0
            
mu = Permeability(markers, degree=1)

# Define current densities
f = Constant(1000.0)/(thickness * width)  #J

# Define variational problem
phi = Function(V)     # rho * A_theta
v = TestFunction(V)

# Poisson equation in axisymmetric cylindrical coordinates
r = Expression('x[0]', degree=1)
F = (1/mu)*(1/r)*(Dx(v,0)*Dx(phi,0) + Dx(v,1)*Dx(phi,1))*dx - f*v*dx(1)

# Compute solution
solve(F==0, phi, bc)
1 Like
2

I’m not really sure what you are trying to do, but you have several parts of your code that is highly inefficient:

calls

which calls a projection over the whole mesh for every cell marked with 2. This does not make sense in any kind of code.
Without having the variational problem you would like to solve in its strong form, it is very hard to give any furher guidance as to how to proceed.

it seems like you want to scale your solution by the previous solution \phi (in the Newton iteration)
f \left(\sqrt{\int_{\Omega} (\frac{1}{x}\frac{\partial \phi}{\partial y})^2+(\frac{1}{x}\frac{\partial \phi}{\partial x} )^2~\mathrm{d}x}\right) (which should only be computed once per iteration).

For such problems, I would strongly suggest setting up your own Picard or Newton solver.
A very old tutorial on how to do this can be found at: A FEniCS Tutorial
and a newer one for Newton’s method using DOLFINx can be found at:
Custom Newton solvers — FEniCSx tutorial

2 Likes
3

Thanks for the help @dokken.

I knew that my code was projecting B over the whole mesh for every cell marked with 2, but I couldn’t figure out another way to do this.
In my problem, there are two subdomains and just one of them is nonlinear (all the cells marked with 2).

Is there any example like this that can help me?

4

As I said, you should build your own Newton solver, then you can create a constant value f_val = f(expr) which starts as 1, and then after the first iteration is equal to the expression above.

You can incorporate this in your variational formulation.

5

Thanks! The following works for me:

# mesh parameters
nmesh = 50 # mesh density
tol = 1e-8 # mesh tolerance 

# geometry parameters
# air box
r_air = 80e-3
z_air = 80e-3

air = Rectangle(Point(tol,-z_air), Point(r_air,z_air))

# Coil
thickness = 25.4e-3
width = 2*25.4e-3
r_inner = 0.5*25.4e-3 #coil inner radius

coil = Rectangle(Point(r_inner,-width/2), Point(r_inner+thickness,width/2))

# Core, the air region inside the coil
core = Rectangle(Point(tol,-width/2), Point(r_inner,width/2))

# domain
domain = air + coil + core

# Set subdomain for wire
domain.set_subdomain(1,coil)
domain.set_subdomain(2,core)

# Create mesh
mesh = generate_mesh(domain, nmesh)

# Steel
Bi_steel = np.array([0, 0.221950, 0.245515, 0.344303, 0.375573, 0.454417, 0.627981, 
                     0.670134, 0.861453, 1.075180, 1.241074, 1.423388, 1.656238, 1.686626, 
                     1.813505, 1.964422, 1.979083, 2.012433, 2.021337, 2.033503, 2.050973, 
                     2.052071, 2.191983, 2.197328, 2.240825, 2.309729, 2.327795, 2.435784])
Hi_steel = np.array([tol, 225.366667, 237.316667, 291.793333, 310.450000, 358.730000, 
                     483.890000, 520.136667, 723.673333, 1071.333333, 1570.566667, 
                     2775.500000, 6290.533333, 7049.866667, 12338.666667, 26304.666667, 
                     28581.000000, 36287.000000, 39022.333333, 43292.333333, 50590.000000, 
                     51118.333333, 134313.333333, 138566.666667, 168803.333333, 223476.666667, 
                     237853.333333, 321480.000000])

mi_steel = np.divide(Bi_steel, Hi_steel)
mi_steel[0] = -Bi_steel[1]*((mi_steel[2] - mi_steel[1])/(Bi_steel[2] - Bi_steel[1])) + mi_steel[1]

f2 = interp1d(Bi_steel, mi_steel, kind='quadratic')

# Define function space
V = FunctionSpace(mesh, 'P', 2)

# Define boundary condition
bc = DirichletBC(V, Constant(0), 'on_boundary')

def mi(phi):
    B = project(as_vector((-(1/r)*phi.dx(1), (1/r)*phi.dx(0))))
    norm_B = norm(B)
    return f2(norm_B)

# Define subdomain markers and integration measure
markers = MeshFunction('size_t', mesh, 2, mesh.domains())
dx = Measure('dx', domain=mesh, subdomain_data=markers)

# Define magnetic permeability
class Permeability(UserExpression):
    def __init__(self, markers, norm_B, **kwargs):
        super().__init__(**kwargs)
        self.markers = markers
        self.norm_B = norm_B
        self.mi_0 = mi0            #coil and air
    
    def eval_cell(self, values, x, cell):
        if self.markers[cell.index] == 2:
            values[0] = f2(self.norm_B([x[0],x[1]]))
        else:
            values[0] = self.mi_0

# Define current densities
f = Constant(1000.0)/(thickness * width)  #J

r = Expression('x[0]', degree=1)
W = VectorFunctionSpace(mesh, 'P', 2)
phi = TrialFunction(V)
v = TestFunction(V)

# Define variational problem for Picard iteration
phi_k = interpolate(Constant(0.0), V)  # previous (known) phi

mu = Permeability(markers,Constant(0.0), degree=1)

a = (1/mu)*(1/r)*(Dx(v,0)*Dx(phi,0) + Dx(v,1)*Dx(phi,1))*dx#*2*pi
L = f*v*dx(1)

# Picard iterations
phi = Function(V)     # new unknown function
eps = 1.0           # error measure ||phi-phi_k||
tol = 1.0E-5        # tolerance
iter = 0            # iteration counter
maxiter = 25        # max no of iterations allowed
while eps > tol and iter < maxiter:
    iter += 1
    solve(a == L, phi, bc)
    dif = phi.vector() - phi_k.vector()
    eps = norm(dif)
    print("iter %d: norm=%g" % (iter, eps))
    phi_k.assign(phi)   # update for next iteration
    B = project(as_vector((-(1/r)*phi.dx(1), (1/r)*phi.dx(0))), W)
    norm_B = sqrt(inner(B,B))
    mu = Permeability(markers,norm_B, degree=1)
6

You are reassigning mu here, and this will not be picked up in a or L.
You should rather do
mu.norm_B = norm_B

1 Like