I want to solve a coupled mechanoc-magnetostatic problem by using a nonlinear solver in a staggered scheme. My domain includes tow copper wires and a rectangular iron. All these are located inside vacuum. The rectangular iron is considered as a single supported beam. I solve the PDE of magnetostatic in whole domain, but I want to solve the mechanical problem for the iron beam. This idea works when I use the linear solver for the mechanical problem, but when I use the nonlinear solver, it does not work and mechanical part is not solved for the beam. Does anybody know solve how to solve this problem?
Regards,
Mehran
from fenics import *
from dolfin import *
from mshr import *
from ufl import nabla_div
import matplotlib.pyplot as plt
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["representation"] = "uflacs"
class Problem(NonlinearProblem):
def __init__(self, J, F, bcs):
self.bilinear_form = J
self.linear_form = F
self.bcs = bcs
NonlinearProblem.__init__(self)
def F(self, b, x):
assemble(self.linear_form, tensor=b)
for bc in self.bcs:
bc.apply(b, x)
def J(self, A, x):
assemble(self.bilinear_form, tensor=A)
for bc in self.bcs:
bc.apply(A)
class CustomSolver(NewtonSolver):
def __init__(self):
NewtonSolver.__init__(self, mesh.mpi_comm(),
PETScKrylovSolver(), PETScFactory.instance())
def solver_setup(self, A, P, problem, iteration):
self.linear_solver().set_operator(A)
PETScOptions.set("ksp_type", "gmres")
# PETScOptions.set("ksp_monitor")
PETScOptions.set("pc_type", "ilu")
self.linear_solver().set_from_options()
T = 10 # final time
num_steps = 2000 # number of time steps
dt = T / num_steps # time step size
alpp = 2.5
alpm = -2.5
pvd_mesh = File("./FEniCSMesh/Mesh.pvd")
### Define the Domain and Subdomains
domain = Rectangle(Point(-0.5,-0.5,), Point(0.5,0.5))
wire_1 = Circle(Point(0.,0.25), 0.1)
domain.set_subdomain(1,wire_1)
rectangle = Rectangle(Point(-0.35, 0.01), Point(0.35,0.11))
domain.set_subdomain(2,rectangle)
wire_2 = Circle(Point(0.,-0.25), 0.1)
domain.set_subdomain(3,wire_2)
### Creat Mesh
mesh = generate_mesh(domain, 50)
mesh_data = open('Mesh.data', 'w')
numel = mesh.num_cells()
mesh_hmax = mesh.hmax()
mesh_hmin = mesh.hmin()
print("Number_of_Elements=", mesh.num_cells())
print("Max_Element_Size=",mesh.hmax())
print("Min_Element_Size=",mesh.hmin())
mesh_data.write( \
"Numer of Elements=%i\n Max. Element Size=%e\n Min. Element Size=%e" \
%(numel, mesh_hmax, mesh_hmin))
mesh_data.close()
pvd_mesh << mesh
# Define function space
V = FunctionSpace(mesh, 'P', 1)
W = VectorFunctionSpace(mesh, 'P', 1)
WW = TensorFunctionSpace(mesh, 'P', 1)
# Define boundary condition
val_init = Constant(0.)
bc = DirichletBC(V, val_init, 'on_boundary')
BCs_emag = [bc]
# Define subdomain markers and integration measure
markers = MeshFunction('size_t', mesh, mesh.topology().dim(), mesh.domains())
plot(markers)
plt.savefig("mesh.png")
dx = Measure('dx', domain=mesh, subdomain_data=markers)
# Define current densities
# J_N = Constant(+1.0)
# J_S = Constant(-1.0)
def curl2D(u):
return as_vector((u.dx(1), -u.dx(0)))
J_N = Expression('alpp*t', degree =1, alpp=alpp, t=0)
J_S = Expression('alpm*t', degree =1, alpm=alpm, t=0)
# Define magnetic permeability
class Permeability(UserExpression):
def __init__(self, markers, **kwargs):
super(Permeability, self).__init__(**kwargs)
self.markers = markers
def eval_cell(self, values, x, cell):
if self.markers[cell.index] == 0:
values[0] = 4*pi*1e-7 # vacuum
elif self.markers[cell.index] == 2:
values[0] = 1e-5 # iron (should really be 6.3e-3)
else:
values[0] = 1.26e-6 # copper
mu = Permeability(markers, degree=1)
#Az_n = TestFunction(V)
Az_n = interpolate(Constant(0.), V)
# Define variational problem
#A_z = TrialFunction(V)
A_z = Function(V)
d = A_z.geometric_dimension() # space dimension
dA_z = TestFunction(V)
a = A_z*dA_z*dx +(1 / mu)*dt*dot(grad(A_z), grad(dA_z))*dx
L_N = J_N*dA_z*dx(1)
L_S = J_S*dA_z*dx(3)
L = L_N + L_S + Az_n*dA_z*dx
F = a - L
J = derivative(F, A_z)
# Solve variational problem
problem = Problem(J, F, BCs_emag)
custom_solver = CustomSolver()
lambda_mech = 2000.
mu_mech = 200.
du = TrialFunction(W)
u = Function(W)
v = TestFunction(W)
dd = u.geometric_dimension()
def epsilon(u):
return 0.5*(nabla_grad(u) + nabla_grad(u).T)
#return sym(nabla_grad(u))
def sigma(u,v):
tau_1 = lambda_mech*nabla_div(u)*Identity(dd) + 2*mu_mech*epsilon(u)
B = curl2D(v)
tau_2 = outer(B,B) + outer(B,B)*epsilon(u)*epsilon(u)
return tau_1 + tau_2
tol_bc = 1E-14
def LeftBC(x, on_boundary):
return near(x[0],-0.35,tol_bc) and near(x[1],+0.01,tol_bc) and on_boundary
def RightBC(x, on_boundary):
return near(x[0],+0.35,tol_bc) and near(x[1],+0.01,tol_bc) and on_boundary
BC_L_x = DirichletBC(W.sub(0), Constant(0.), LeftBC, method ='pointwise')
BC_L_y = DirichletBC(W.sub(1), Constant(0.), LeftBC, method ='pointwise')
BC_R_y = DirichletBC(W.sub(1), Constant(0.),RightBC, method ='pointwise')
BCs_mech = [BC_L_x, BC_L_y, BC_R_y]
F_mech = inner(sigma(u, Az_n), epsilon(v))*dx(2)
a_mech = lhs(F_mech)
L_mech = rhs(F_mech)
J_mech = derivative(F_mech, u)
problem_mech = Problem(J_mech, F_mech, BCs_mech)
custom_solver_mech = CustomSolver()
vtkfile_A_z = File('./Results/A_z.pvd')
vtkfile_B = File('./Results/B.pvd')
vtkfile_disp = File('./Results/Displacement.pvd')
t=0.
print("dTime=", dt)
for n in range(num_steps):
t+=dt
J_N.t = t
J_S.t = t
custom_solver.solve(problem, A_z.vector())
custom_solver_mech.solve(problem_mech, u.vector())
# Compute magnetic field (B = curl A)
B = curl2D(A_z)
B_Field = project(B, W)
# Save solution to file
A_z.rename("A_z","1")
B_Field.rename("B","2")
u.rename("Displacement","3")
vtkfile_A_z << A_z
vtkfile_B << B_Field
vtkfile_disp << u
Az_n.assign(A_z)