### Solver Problem in Nonlinear Magnetostatics

291

views

1

Hi everyone,

I'm currently working on nonlinear magnetostatics. I'm using a convenient scalar potential formulation.

I want to use a constitutive model which takes into account the effect of the magnetic saturation, but I'm facing problems trying to solve the nonlinear variational problem.

Here is a short working example to reproduce the error:

```
from dolfin import *
# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["representation"] = "uflacs"
parameters["form_compiler"]["quadrature_degree"] = 4
# Create mesh and define function space
mesh = UnitSquareMesh(32, 32)
V = FunctionSpace(mesh, "Lagrange", 1)
# Mark boundary subdomians
left = CompiledSubDomain("near(x[0], side) && on_boundary", side = 0.0)
right = CompiledSubDomain("near(x[0], side) && on_boundary", side = 1.0)
# Dirichlet boundary::
bcl = DirichletBC(V, Constant(0), left)
bcr = DirichletBC(V, Constant(100e+3), right)
bcs = [bcl, bcr]
# Define functions
du = TrialFunction(V)
v = TestFunction(V)
u = Function(V)
# Kinematics
H = -grad(u)
# Invariant
IN = dot(H, H)
# parameters
mu0= 4.0*pi * 1e-13
phi=10
ms= 1e+6
a=mu0*ms*ms/phi
b=phi/ms
# Magnetic induction
B=mu0*H + a*b*tanh(b*(IN)**0.5)/((IN)**0.5)*H
# Variational Formulation
F = dot(B, grad(v))*dx
# Compute Jacobian of F
J = derivative(F, u, du)
# Solve variational problem
problem = NonlinearVariationalProblem(F, u, bcs, J)
solver = NonlinearVariationalSolver(problem)
solver.parameters['newton_solver']['linear_solver'] = 'mumps'
solver.solve()
# Save solution in VTK format
file = File("potential.pvd");
file << u;
```

It seems like there is a linearization problem with the expression accounting for the saturation of the magnetization.

I will be grateful for any help.

Best regards,

Dustin

Community: FEniCS Project

### 2 Answers

4

Hi there,

please be cautious about dividing by the invariant, which is zero (initially). Use something like: IN=dot(H,H)+1E-10 in order to circumvent this issue.

Best, Emek

please be cautious about dividing by the invariant, which is zero (initially). Use something like: IN=dot(H,H)+1E-10 in order to circumvent this issue.

Best, Emek

0

Hi Emek,

thanks for the little hint. Of course this was quite obvious ...

Best regards,

Dustin

thanks for the little hint. Of course this was quite obvious ...

Best regards,

Dustin

Please login to add an answer/comment or follow this question.

`u`

?