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 form the mechanical problem, but when I use the nonlinear solver, it does not work and any displacement solved for the beam. Does anybody know solve this problem?
Best Regards,
Mehran
The code is below:
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)