Dear colleagues,
I want to solve Usadel equations with fenicsx for a structure superconductor/ ferromagnet/ superconductor. These equations are formulated for 8 complex functions, \chi^{1}, … \chi^{8} with the help of 2to2 matrices:
\hat{\gamma} = \chi^1 \hat{1} + \chi^{2} \hat{\sigma}_x + \chi^{3} \hat{\sigma}_y + \chi^{4} \hat{\sigma}_z,\qquad
\hat{ \tilde{\gamma} } = \chi^5 \hat{1} + \chi^{6} \hat{\sigma}_x + \chi^{7} \hat{\sigma}_y + \chi^{8}\hat{\sigma}_z;\\
\hat{N} = (\hat{1} + \hat{\gamma}\hat{\tilde{\gamma}})^{-1}, \qquad \hat{\tilde{N}} = -(\hat{1} + \hat{\tilde{\gamma}}\hat{\gamma})^{-1}.\\
Here hats \hat{\ } are used for matrices 2 to 2, \hat{\sigma}_{x,y,z} are Pauli matrices.
Using there matrices, I formulate Usadel equations as \nabla^2 \hat{\gamma} + 2\nabla\hat{\gamma} \hat{\tilde{N}}\hat{\tilde{\gamma}} \nabla{\hat{\gamma}} = \omega_n \hat{\gamma},\ \nabla^2 \hat{\tilde{\gamma}} - 2\nabla\hat{\tilde{\gamma}} \hat{{N}}\hat{\gamma} \nabla{\hat{\tilde{\gamma}}} = \omega_n \hat{\tilde{\gamma}}.
These equations are non-linear, to solve them I use Picard iterations which go as follows:
\hat{\tilde{N}}^{(n)} = -\left(1 + \hat{\tilde{\gamma}}^{(n)} \hat{\tilde{\gamma}}^{(n)}\right)^{-1}, \qquad \hat{N}^{(n)} = \left(1 + \hat{\gamma}^{(n)}\hat{\tilde{\gamma}}^{(n)}\right)^{-1};\\
2 \nabla' \hat{\gamma}\hat{\tilde{N}} \hat{\tilde{\gamma}} \nabla\hat{\gamma} \to \frac{1}{3}\left\{ 2 \nabla' \hat{\gamma}^{(n+1)}\hat{\tilde{N}}^{(n)} \hat{\tilde{\gamma}}^{(n)} \nabla\hat{\gamma}^{(n)} + 2 \nabla' \hat{\gamma}^{(n)}\hat{\tilde{N}}^{(n)} \hat{\tilde{\gamma}}^{(n)} \nabla\hat{\gamma}^{(n+1)} + 2 \nabla' \hat{\gamma}^{(n)}\hat{\tilde{N}}^{(n)} \hat{\tilde{\gamma}}^{(n+1)} \nabla\hat{\gamma}^{(n)} \right\},\\
2 \nabla' \hat{\tilde{\gamma}}\hat{N} \hat{\gamma} \nabla\hat{\tilde{\gamma}} \to \frac{1}{3}\left\{ 2 \nabla' \hat{\tilde{\gamma}}^{(n+1)}\hat{N}^{(n)} \hat{\gamma}^{(n)} \nabla\hat{\tilde{\gamma}}^{(n)} + 2 \nabla' \hat{\tilde{\gamma}}^{(n)}\hat{N}^{(n)} \hat{\gamma}^{(n)} \nabla\hat{\tilde{\gamma}}^{(n+1)} + 2 \nabla' \hat{\tilde{\gamma}}^{(n)}\hat{N}^{(n)} \hat{\gamma}^{(n+1)} \nabla\hat{\tilde{\gamma}}^{(n)} \right\}.
Here n is the number of the current iteration in non-linearity. After I obtain the solution of the full Usadel equations after a few iterations in the ferromagnet, I calculate charge fluxes, one from the left superconductor to the ferromagnet, I_L, the another from the ferromagnet to superconductor, I_R. Superconductors play a role of boundary conditions only.
These currents are calculated as integrals: I \propto \int dS \operatorname{Im}\ \operatorname{Tr}_2\left\{ \left(1 + \hat{\tilde{\gamma}}\hat{\gamma}\right)^{-3}\hat{\tilde{\gamma}}\left[\nabla\hat{{\gamma}} - \hat{{\gamma}}(\nabla \hat{\tilde{\gamma}})\hat{{\gamma}}\right] \right\}
From physical point of view, I have to obtain the relation I_L = I_R. However, for a MixedElement containing linear functions, I didn’t succeeded. As far as I understand, the problem comes from too small degree of Lagrange polynomials. However, the program appears to be too slow, hours. Are there any ideas how to speed the code up?
MWE goes as follows:
from mpi4py import MPI
import ufl
import dolfinx
import dolfinx.fem.petsc as petsc
import numpy as np
import sympy as sp
from sympy.physics import msigma
import time
t_start = time.time()
# dimensionless geometrical sizes
d_F = 20./25.
L_y = 100./25.
w_F = 30./25.
d_S = 50./25.
t = 0.1 # temperature t = T/T_c
N_F_x = 10*int(w_F)
N_F_y = 10*int(L_y)
N_F_z = 10*int(d_F +2*d_S)
print(" The mesh with the dimensions N_F = ", (N_F_x, N_F_y, N_F_z)," is beign created...")
mesh_In_F = dolfinx.mesh.create_box(comm=MPI.COMM_WORLD,
points=((-w_F, L_y/2, -d_S - d_F/2), (0., -L_y/2, d_F/2 + d_S)),
n=(N_F_x, N_F_y, N_F_z),
cell_type=dolfinx.mesh.CellType.tetrahedron)
chi_el = ufl.FiniteElement("CG", mesh_In_F.ufl_cell(), 1)
V = dolfinx.fem.FunctionSpace(mesh_In_F, ufl.MixedElement([chi_el,
chi_el,
chi_el,
chi_el,
chi_el,
chi_el,
chi_el,
chi_el]))
print(" Mesh was created. ")
def is_left_lead(p):
return np.logical_and.reduce((-d_S - d_F/2 <= p[2], p[2] <= - d_F/2, - L_y/2 <= p[1], p[1] <= L_y/2, np.isclose(p[0], 0.) ))
def is_right_lead(p):
return np.logical_and.reduce((d_F/2 <= p[2], p[2] <= d_F/2 + d_S, - L_y/2 <= p[1], p[1] <= L_y/2, np.isclose(p[0], 0.) ))
mesh_In_F.topology.create_connectivity(mesh_In_F.topology.dim-1, mesh_In_F.topology.dim)
SF_R_facets = dolfinx.mesh.locate_entities_boundary(mesh_In_F, mesh_In_F.topology.dim - 1, is_right_lead)
SF_R_facets_mark = np.zeros_like(SF_R_facets) + 1
SF_R_subdomain_data = dolfinx.mesh.meshtags(mesh_In_F, mesh_In_F.topology.dim-1, np.array(
SF_R_facets).astype(np.int32), np.array(SF_R_facets_mark).astype(np.int32))
SF_L_facets = dolfinx.mesh.locate_entities_boundary(mesh_In_F, mesh_In_F.topology.dim - 1, is_left_lead)
SF_L_facets_mark = np.zeros_like(SF_L_facets) + 2
SF_L_subdomain_data = dolfinx.mesh.meshtags(mesh_In_F, mesh_In_F.topology.dim-1, np.array(
SF_L_facets).astype(np.int32), np.array(SF_L_facets_mark).astype(np.int32))
dofs_L = [dolfinx.fem.locate_dofs_geometrical((V.sub(i), V.sub(i).collapse()[0]), is_left_lead) for i in range(V.num_sub_spaces)]
dofs_R = [dolfinx.fem.locate_dofs_geometrical((V.sub(i), V.sub(i).collapse()[0]), is_right_lead) for i in range(V.num_sub_spaces)]
# here iter is the index of the current iteration in non-linearity for the Usadel equation
def getBCs(chi_L, chi_R):
bc_L = [dolfinx.fem.dirichletbc(dolfinx.fem.Constant(mesh_In_F, petsc.PETSc.ScalarType(chi_L[i])), dofs_L[i][0], V.sub(i)) for i in range(V.num_sub_spaces)]
bc_R = [dolfinx.fem.dirichletbc(dolfinx.fem.Constant(mesh_In_F, petsc.PETSc.ScalarType(chi_R[i])), dofs_R[i][0], V.sub(i)) for i in range(V.num_sub_spaces)]
BCs = [*bc_L, *bc_R]
return BCs
omega_n = dolfinx.fem.Constant(mesh_In_F, petsc.PETSc.ScalarType((2*0 + 1)*t))
Delta_s = 1.
varphi_L = -np.pi/3
varphi_R = np.pi/3
chi_L = [0, 0, np.exp(1j*varphi_L)*(omega_n.value/Delta_s - np.sqrt(omega_n.value**2/Delta_s**2 + 1.)), 0,
0, 0, np.exp(-1j*varphi_L)*(omega_n.value/Delta_s - np.sqrt(omega_n.value**2/Delta_s**2 + 1.)), 0]
chi_R = [0, 0, np.exp(1j*varphi_R)*(omega_n.value/Delta_s - np.sqrt(omega_n.value**2/Delta_s**2 + 1.)), 0,
0, 0, np.exp(-1j*varphi_R)*(omega_n.value/Delta_s - np.sqrt(omega_n.value**2/Delta_s**2 + 1.)), 0]
BCs = getBCs(chi_L, chi_R)
eta = ufl.TestFunction(V)
chi_TF = ufl.TrialFunction(V)
chi_prev = dolfinx.fem.Function(V, dtype=np.complex128)
for i in range(V.num_sub_spaces):
chi_prev.sub(i).interpolate(lambda x: np.full((x.shape[1],), 0.))
chi = dolfinx.fem.Function(V, dtype=np.complex128)
for i in range(V.num_sub_spaces):
chi.sub(i).interpolate(lambda x: np.full((x.shape[1],), 1. ))
sigma_basis_p = [sp.eye(2), msigma(1), msigma(2), msigma(3)]
sigma_basis = [ufl.as_matrix(np.array(sigma_basis_p[i], dtype=np.complex128).tolist()) for i in range(len(sigma_basis_p)) ]
gamma_m = ufl.as_matrix(sigma_basis[0]*chi_prev[0] + sigma_basis[1]*chi_prev[1] +
sigma_basis[2]*chi_prev[2] + sigma_basis[3]*chi_prev[3])
gamma_mt = ufl.as_matrix(sigma_basis[0]*chi_prev[4] + sigma_basis[1]*chi_prev[5] +
sigma_basis[2]*chi_prev[6] + sigma_basis[3]*chi_prev[7])
gamma_mx = ufl.as_matrix(sigma_basis[0]*chi_TF[0] + sigma_basis[1]*chi_TF[1] +
sigma_basis[2]*chi_TF[2] + sigma_basis[3]*chi_TF[3])
gamma_mtx = ufl.as_matrix(sigma_basis[0]*chi_TF[4] + sigma_basis[1]*chi_TF[5] +
sigma_basis[2]*chi_TF[6] + sigma_basis[3]*chi_TF[7])
N_m = ufl.inv(ufl.Identity(2) + gamma_m*gamma_mt)
N_mt = - ufl.inv(ufl.Identity(2) + gamma_mt*gamma_m)
A = +2*( ufl.Dx(gamma_mx, 0)*N_mt*gamma_mt*ufl.Dx(gamma_m, 0) +\
ufl.Dx(gamma_mx, 1)*N_mt*gamma_mt*ufl.Dx(gamma_m, 1) +\
ufl.Dx(gamma_mx, 2)*N_mt*gamma_mt*ufl.Dx(gamma_m, 2) +\
ufl.Dx(gamma_m, 0)*N_mt*gamma_mt*ufl.Dx(gamma_mx, 0) +\
ufl.Dx(gamma_m, 1)*N_mt*gamma_mt*ufl.Dx(gamma_mx, 1) +\
ufl.Dx(gamma_m, 2)*N_mt*gamma_mt*ufl.Dx(gamma_mx, 2) +\
ufl.Dx(gamma_m, 0)*N_mt*gamma_mtx*ufl.Dx(gamma_m, 0) +\
ufl.Dx(gamma_m, 1)*N_mt*gamma_mtx*ufl.Dx(gamma_m, 1) +\
ufl.Dx(gamma_m, 2)*N_mt*gamma_mtx*ufl.Dx(gamma_m, 2) )/3
B = -2*( ufl.Dx(gamma_mtx, 0)*N_m*gamma_m*ufl.Dx(gamma_mt, 0) +\
ufl.Dx(gamma_mtx, 1)*N_m*gamma_m*ufl.Dx(gamma_mt, 1) +\
ufl.Dx(gamma_mtx, 2)*N_m*gamma_m*ufl.Dx(gamma_mt, 2) +\
ufl.Dx(gamma_mt, 0)*N_m*gamma_m*ufl.Dx(gamma_mtx, 0) +\
ufl.Dx(gamma_mt, 1)*N_m*gamma_m*ufl.Dx(gamma_mtx, 1) +\
ufl.Dx(gamma_mt, 2)*N_m*gamma_m*ufl.Dx(gamma_mtx, 2) +\
ufl.Dx(gamma_mt, 0)*N_m*gamma_mx*ufl.Dx(gamma_mt, 0) +\
ufl.Dx(gamma_mt, 1)*N_m*gamma_mx*ufl.Dx(gamma_mt, 1) +\
ufl.Dx(gamma_mt, 2)*N_m*gamma_mx*ufl.Dx(gamma_mt, 2) )/3
def proj_2_2_to_8_1(A, B):
sigma_basis_p = [sp.eye(2), msigma(1), msigma(2), msigma(3)]
sigma_basis = [ufl.as_matrix(np.array(sigma_basis_p[i], dtype=np.complex128).tolist()) for i in range(len(sigma_basis_p)) ]
return [ufl.tr(A*sigma_basis[0])/2,
ufl.tr(A*sigma_basis[1])/2,
ufl.tr(A*sigma_basis[2])/2,
ufl.tr(A*sigma_basis[3])/2,
ufl.tr(B*sigma_basis[0])/2,
ufl.tr(B*sigma_basis[1])/2,
ufl.tr(B*sigma_basis[2])/2,
ufl.tr(B*sigma_basis[3])/2]
def get_norm_of_difference(a, b):
comm = a.function_space.mesh.comm
difference = dolfinx.fem.form((a - b)**2 * ufl.dx)
E = np.abs(np.sqrt(comm.allreduce(dolfinx.fem.assemble_scalar(difference), MPI.SUM)))
return E
def F0_mm_f_ufl(x, omega_n):
return omega_n*ufl.Identity(8)
x = ufl.SpatialCoordinate(mesh_In_F)
LeftHS = - ufl.inner( ufl.grad(chi_TF), ufl.grad(eta)) + \
ufl.inner(ufl.as_vector(proj_2_2_to_8_1(A, B)), eta) - \
ufl.inner(F0_mm_f_ufl(x, omega_n)*chi_TF, eta)
dx = ufl.Measure("dx", domain=mesh_In_F)
LeftHS_fem_form = dolfinx.fem.form(LeftHS*dx)
I_n = dolfinx.fem.Constant(mesh_In_F, petsc.PETSc.ScalarType((0, 0, 0, 0, 0, 0, 0, 0)))
RightHS_fem_form = dolfinx.fem.form(ufl.inner(I_n, eta)*dx)
omega_n.value = (2*0 + 1)*t
LeftHS_mmm = dolfinx.fem.petsc.assemble_matrix(LeftHS_fem_form, bcs=BCs)
LeftHS_mmm.assemble()
RightHS_v = dolfinx.fem.petsc.create_vector(RightHS_fem_form)
solver = petsc.PETSc.KSP().create(mesh_In_F.comm)
solver.setType(petsc.PETSc.KSP.Type.GMRES)
solver.getPC().setType(petsc.PETSc.PC.Type.ILU)
solver.setOperators(LeftHS_mmm)
iter = 0
norm_of_difference = get_norm_of_difference(chi, chi_prev)
iter_max = 25
atol = 1e-4
print(" Iterations in non-linearity: ")
while (iter <= iter_max and norm_of_difference >= atol):
with RightHS_v.localForm() as loc_1:
loc_1.set(0)
RightHS_v= dolfinx.fem.petsc.assemble_vector(RightHS_v, RightHS_fem_form)
dolfinx.fem.petsc.apply_lifting( RightHS_v, [LeftHS_fem_form], [BCs])
RightHS_v.ghostUpdate(addv=petsc.PETSc.InsertMode.ADD_VALUES, mode=petsc.PETSc.ScatterMode.REVERSE)
dolfinx.fem.petsc.set_bc(RightHS_v, BCs)
solver.solve(RightHS_v, chi.vector)
chi.x.scatter_forward()
norm_of_difference = get_norm_of_difference(chi, chi_prev)
chi_prev.x.array[:] = chi.x.array
chi_prev.x.scatter_forward()
iter = iter + 1
LeftHS_mmm = dolfinx.fem.petsc.assemble_matrix(LeftHS_fem_form, bcs=BCs)
LeftHS_mmm.assemble()
solver.setOperators(LeftHS_mmm)
print(" Iter = ", iter, " is finishesd with norm_of_difference =", norm_of_difference)
gamma_m = ufl.as_matrix(sigma_basis[0]*chi[0] + sigma_basis[1]*chi[1] +
sigma_basis[2]*chi[2] + sigma_basis[3]*chi[3])
dgamma_mx = ufl.as_matrix(sigma_basis[0]*ufl.Dx(chi[0], 0) + sigma_basis[1]*ufl.Dx(chi[1], 0) +
sigma_basis[2]*ufl.Dx(chi[2], 0) + sigma_basis[3]*ufl.Dx(chi[3], 0))
gamma_mt = ufl.as_matrix(sigma_basis[0]*chi[4] + sigma_basis[1]*chi[5] +
sigma_basis[2]*chi[6] + sigma_basis[3]*chi[7])
dgamma_mtx = ufl.as_matrix(sigma_basis[0]*ufl.Dx(chi[4], 0) + sigma_basis[1]*ufl.Dx(chi[5], 0) +
sigma_basis[2]*ufl.Dx(chi[6], 0) + sigma_basis[3]*ufl.Dx(chi[7], 0))
N_mt = - ufl.inv(ufl.Identity(2) + gamma_mt*gamma_m)
ds = ufl.Measure("ds", domain=mesh_In_F, subdomain_data=SF_R_subdomain_data)
j_x = dolfinx.fem.form( 4*ufl.pi*ufl.tr(-N_mt*N_mt*N_mt*gamma_mt*( dgamma_mx - gamma_m*dgamma_mtx*gamma_m) )*ds(1) )
comm = chi.function_space.mesh.comm
I = t*np.imag(comm.allreduce(dolfinx.fem.assemble_scalar(j_x), MPI.SUM))
print(" currenr from F to right S: ", I)
ds = ufl.Measure("ds", domain=mesh_In_F, subdomain_data=SF_L_subdomain_data)
j_x = dolfinx.fem.form( 4*ufl.pi*ufl.tr(-N_mt*N_mt*N_mt*gamma_mt*( dgamma_mx - gamma_m*dgamma_mtx*gamma_m) )*ds(2) )
comm = chi.function_space.mesh.comm
I = t*np.imag(comm.allreduce(dolfinx.fem.assemble_scalar(j_x), MPI.SUM))
print(" currenr from left S to F: ", I)
print(" time used = ", time.time() - t_start)
