### Solving the Ginzburg-Landau-Equations on UnitSquareMesh and quadratic domain gives different results

292
views
1
10 months ago by
Dear FEniCS Community,

I have the rather difficult task to solve the time dependent Ginzburg-Landau-Equations
$\frac{\partial \psi}{\partial t} - i \kappa ( \nabla \cdot \vec A) \psi +\left(\frac{i}{\kappa} \nabla + \vec A \right)^2 \psi +\left(|\psi|^2-1 \right)\psi =0 \\ \frac{\partial \vec A}{\partial t}-\nabla (\nabla \cdot \vec A) +\nabla \times \vec \sigma +\frac{i}{2 \kappa} \left(\psi^{*} \nabla \psi- \psi \nabla \psi^{*} \right) +|\psi|^2 \vec A = \nabla \times \vec H \\ \sigma = \nabla \times \vec A$

with given boundary conditions and initial conditions
$\vec A \cdot \vec n =0 ~~~~~~ \text{on \Gamma \times [0,T]} \\ \psi =0 ~~~~~~ \text{on \Gamma \times [0,T]} \\ \vec \sigma \times \vec n=\vec H \times \vec n ~~~~~~ \text{on \Gamma \times [0,T]} \\ \psi(x,0)=\psi_{0}(x) ~~~~~~ \text{in \Omega} \\ \vec A (x,0)=\vec A_{0} (x) ~~~~~~ \text{in \Omega}$
A mixed finite element methode is used to solve the equations in two dimensions for $\psi$ and $( \sigma ,\vec A)$ where $\psi =(\psi_{re}, \psi_{im})$ is a complex function, $\sigma$ and $H$ are scalars and $\vec A$ is a vector field.  The domain of computation is the unit square. Now the problem is that it is possible to choose a symmetric partition with UnitSquareMesh(n,n) or a nonsymmetric one by defining a domain with Rectangle(Point(0,0),Point(1,1)). The solution for $|\psi|$ on UnitSquareMesh however is less symmetric than the solution on the other domain. this should not be the case. I provide a minimal running example to demonstrate what I mean. Another problem with the partitions is that the function $\Phi = - \nabla \cdot \vec A$ is less smooth on the rectangle-domain and becomes really rough for more complicated fields  $H$ . Has someone an idea why there are different solutions on the same domain while all other parameters are the same? Thanks!

from dolfin import *
import numpy as np
import scipy
import math
from mshr import  *

#constants and H(x[0], x[1], tc)==================================================

T = 50        #final time
TOL = 3E-15   #parameter for the solver
dt = 1.0      #width of the time steps
tc = dt       #current time
kappa = 100   #parameter

H = Expression('10', degree=3)# external magnetic field

#mesh=============================================================================
mesh = UnitSquareMesh(80,80)

#domain = Rectangle(Point(0,0),Point(1,1))
#mesh = generate_mesh(domain,80)

#Function spaces and trial/test functions==========================================
V1 = FunctionSpace(mesh,'CG',2)#          Function space for the magnetic field H
W1 = VectorFunctionSpace(mesh, 'CG', 1)#  Function space for curl(H)
U = VectorFunctionSpace(mesh, "CG", 1)#   Function space for the order parameter psi
W =  FunctionSpace(mesh, "RT", 1)#        Function space for the Vector potential A
VA = FunctionSpace(mesh, "CG", 1)#        Function space for sigma = curl(A)

RT = FiniteElement("RT", mesh.ufl_cell(), 1)   #element for A
CG  = FiniteElement("CG", mesh.ufl_cell(), 1)  #element for sigma
M = FunctionSpace(mesh, CG * RT)#               Mixed Functionspace for the vector (sigma, A)

(psi_0,psi_1) = TrialFunctions(U)
(v0,v1) = TestFunctions(U)

(sigma,A) = TrialFunctions(M)
(tau,v) = TestFunctions(M)

BIG1 = Function(M)
sigma1 = Function(VA)

#initial functions and boundary values=========================================

def boundary(x,on_boundary):
return on_boundary

k = Expression("{}".format(kappa),degree=3)

p0_e = Expression(("0.8","0.6"),degree=3)# Initial funcction for psi
A0_e = Expression(("0","0"),degree=3)    # Initial function for A

psi_b = Constant((0,0))                  # boundary value of psi

p0 = interpolate(p0_e,U)                 # psi in time step n-1
p1 = interpolate(p0_e,U)                 # psi in time step n
(ph0_0, ph0_1) = p0.split()
(ph1_0, ph1_1) = p1.split()

Ah0 = Function(W)
Ah1 = Function(W)

#boundary condition for A-----------------------------------------------------
class BoundarySource(Expression):
def __init__(self, mesh, **kwargs):
self.mesh = mesh
def eval_cell(self, values, x, ufc_cell):
cell = Cell(self.mesh, ufc_cell.index)
n = cell.normal(ufc_cell.local_facet)
g = 0
values[0] = g*n[0]
values[1] = g*n[1]
def value_shape(self):
return (2,)

G = BoundarySource(mesh, degree=2)
#-------------------------------------------------------------------------

#magnetic field========================================================
Hz = project(H, V1)
H_c = project(as_vector((Hz.dx(1),-Hz.dx(0))), W1) #curl(Hz)

#loop==================================================================
while (tc < T + 1e-10):

#===========================
#= solve psi
#===========================

a1 = (psi_0*v0 + psi_1*v1)*dx + \
- dt*k*div(Ah0)*(psi_0*v1 - psi_1*v0)*dx \
+ dt*(dot(Ah0,Ah0) + dot(p0,p0) -1.0)*(psi_0*v0 + psi_1*v1)*dx

L1 = (ph0_0*v0 + ph0_1*v1)*dx

#----------
solverp = KrylovSolver("bicgstab", "amg")
solverp.parameters["absolute_tolerance"] = TOL
solverp.parameters["relative_tolerance"] = TOL
solverp.parameters["maximum_iterations"] = 10000000
solverp.parameters["nonzero_initial_guess"] = True
#----------

bc = DirichletBC(M.sub(0), Hz, boundary)   #boundary condition for sigma

bcAn = DirichletBC(M.sub(1), G, boundary)  #boundary condition for A

bc_psi=DirichletBC(U,psi_b,boundary)       #boundary condition for psi

A_psi = assemble(a1)
b_psi = assemble(L1)

bc_psi.apply(A_psi,b_psi)

solverp.solve(A_psi, p1.vector(), b_psi)

#===========================
#= solve A
#===========================

a2 = (dot(sigma,tau)-dot(curl(tau), A)-dot(curl(sigma),v)-div(A)*div(v))*dx \
- 1.0/dt*dot(A,v)*dx - dot(A,v)*dot(p0,p0)*dx

L2 = - 1.0/dt*dot(Ah0,v)*dx \

A_A = assemble(a2)
b_A = assemble(L2)

bc.apply(A_A,b_A)
bcAn.apply(A_A,b_A)

solverA = LinearSolver("lu")
solverA.solve(A_A, BIG1.vector(), b_A)
(sigma1, Ah1) = BIG1.split(True)

assign(sigma1, BIG1.sub(0))
assign(    Ah1, BIG1.sub(1))
assign(    p0,          p1)
assign(    Ah0,          Ah1)

assign(ph1_0, p1.sub(0))
assign(ph1_1, p1.sub(1))
assign(ph0_0, p0.sub(0))
assign(ph0_1, p0.sub(1))

print( "TC = %12.5f" % tc)

#updating functions
tc = tc + dt
H.tc = tc
Hz = project(H, V1)
H_c = project(as_vector((Hz.dx(1),-Hz.dx(0))), W1)

Vphi=FunctionSpace(mesh,'CG', 1)            # function space for phi
phi=project(Ah1[0].dx(0)+Ah1[1].dx(1),Vphi) # phi

plot(mesh)
plot(dot(p1,p1),mesh,title='psi')
plot(phi,mesh,title='phi')

interactive()​
Community: FEniCS Project