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


292
views
1
10 months ago by
BN319  
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*(1.0/k/k)*(inner(grad(psi_0),grad(v0))+inner(grad(psi_1),grad(v1)))*dx \
        - dt*k*div(Ah0)*(psi_0*v1 - psi_1*v0)*dx \
        + dt*(1.0/k)*(inner(grad(psi_0), Ah0)*v1 - inner(grad(psi_1), Ah0)*v0)*dx \
        - dt*(1.0/k)*(inner(grad(v1), Ah0)*psi_0 - inner(grad(v0), Ah0)*psi_1)*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 \
        + (1.0/k)*(ph0_1*dot(grad(ph0_0),v)-ph0_0*dot(grad(ph0_1),v))*dx-dot(H_c,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
Please login to add an answer/comment or follow this question.

Similar posts:
Search »