Time-dependent Ginzburg-Landau equation

1

Hey Guys,
I hope you are great. I am completely new to Fenics and this is my first time posting to the group. I appreciate that if you can help me.
I am working on some phase field problem. For start, I need to solve a time-dependent Ginzburg-Landau equation with a highly nonlinear source term:
(1/lmda) du/dt = lmbdabettagrad^2(u) - 2Almbda*(u+2u3-3u2)+12lmbdadelG*(u2 - u3)
where all the coefficients are given below in the code. In fact, the RHS is the driving force for propagation of interface. u is in fact a scalar parameter. The initial condition for u is pow(1+exp(x[0] - 0.5), -1).
I have started from Cahn-Hilliard equation, then try the nonlinear time-dependent burgers equation, However I was not able to solve the problem. The latest code that I have written is given below for your consideration:

from fenics import *

Class for interfacing with the Newton solver

class GinzburgLandau(NonlinearProblem):
def init(self, a, L):
NonlinearProblem.init(self)
self.L = L
self.a = a
def F(self, b, x):
assemble(self.L, tensor=b)
def J(self, A, x):
assemble(self.a, tensor=A)

Create mesh and build function space

mesh = UnitSquareMesh(80, 80)
V = FunctionSpace(mesh, ‘P’, 2)

Define trial and test functions

du = TrialFunction(V)
v = TestFunction(V)

Define functions and initial condition

u = Function(V) # current solution
u_n = Function(V) # solution from previous converged step

u_init = Expression(‘pow((1+exp((x[0] - 0.5) / 2.95e-9)), -1)’, degree=2)
u.interpolate(u_init)
u_n.interpolate(u_init)

Weak statement of the equations

lmbda = 2600.0
betta = 2.59e-10
delG = -1073.3e6
A = 7924.0e6
dt = 0.1e-14
k = Constant(dt)

L0 = uvdx - u_nvdx + lmbdabettakdot(grad(u), grad(v))dx
L1 = (2kAlmbda(u + 2u**3 - 3u2) + 12klmbdadelG(u2 - u**3))vdx
L = L0 + L1

Compute directional derivative about u in the direction of du (Jacobian)

a = derivative(L, u, du)

Create nonlinear problem and Newton solver

problem = GinzburgLandau(a, L)
solver = NewtonSolver()
solver.parameters[“linear_solver”] = “lu”

solver.parameters[“preconditioner”] = “ilu”

solver.parameters[“convergence_criterion”] = “incremental”
solver.parameters[“relative_tolerance”] = 1e-12

solver.parameters[“absolute_tolerance”] = 1e-13

solver.parameters[“maximum_iterations”] = 25

solver.parameters[“relaxation_parameter”] = 1.2

file_u = File(‘ginzburglandau/u.pvd’)

num_steps = 500
t = 0.0
for n in range(num_steps):

# Update current time
t += dt
print('t = ', t)
solver.solve(problem, u.vector())
file_u << u

# Update the previous solution
u_n.assign(u)

Unfortunately, the interface does not propagate with time and it stops in the middle of the square. I appreciate that if you can show me a start point that I can get to the solution of my problem.
I know you are really busy, but it would be very kind of you if you can respond it.
Looking forward to hearing form you.

Bests,
Benhour

1 Like
2

Looking at the values of the parameters, I think you’re probably just under-resolving the problem in space (by many orders of magnitude). Doing some dimensional analysis, you can see that \beta / A will have units of length squared, so, heuristically, your mesh element diameter would need to be comparable to \sqrt{\beta / A} to resolve the physics. If I manually pump up \beta to 8\times 10^6, then the interface moves. (Of course this changes the physics of the problem; hopefully you’re instead planning to solve the PDE on a much smaller domain, which would have a similar effect.)

P.S. You can format code by surrounding it with backticks (and optionally specifying the language for highlighting), e.g.,

```python
Line 1
Line 2
etc.
```

For more formatting tricks, Google for “Markdown”. You can also embed equations with LaTeX notation, e.g., $\beta / A$.

3

Dear David,
Thanks very much for your response. Lets talk about a time-dependent nonlinear problem. The problem that I want to solve is:
\frac{\partial u}{\partial t} = \Delta{u} - u - 2u^{3} + 3u^{2} (u is a scalar variable)
with the initial condition for u = \frac{1 + exp(\frac{x - 0.5}{2.95e-9})}{-1}, and homogeneous Neumann boundary condition. By solving this problem, I should have an interface right in the middle which propagates to the left or right (depend on the sign of RHS). Let say the geometry is a unit square. I wrote the following code:

from fenics import *
from mshr import *
%matplotlib inline

R = Rectangle (Point(0, 0), Point(1.0, 1.0))
mesh = generate_mesh(R, 100)
plot(mesh)

V = FunctionSpace(mesh, 'P', 2)
u_init = Expression('pow(1 + exp((x[0] - 0.5) / 2.95e-9), -1)', degree = 2)
t = 0.0
T = 0.1
ts = 0.0001

u = Function(V)
u_n = interpolate(u_init, V)
du = TrialFunction(V)
v = TestFunction(V)
F = (((u - u_n)/ts)*v*dx + inner(grad(u), grad(v))*dx + (u + 2*u*u*u - 3*u*u)
J = derivative(F, u, du)

file1 = File('results/u.pvd')
while t < T
      solve(F == 0, u, J=J)
      file1 << (u, t)
      print('t = ', t)
      u_n.assign(u)
      t+=ts

When I solved this problem, the interface propagates from the center in both direction (decays instead of propagate in only one direction). I have attached some of my results for your consideration. Please let me know where I am making mistakes. In fact the interface profile should be moved to the left or right, not to expand. I have also made my geometry small and then played with the time steps, however I could not get the desired result. I appreciate your time and help if you can take a look at the code and see where it is going wrong.

Bests,
Benhour


![1|410x500](upload://b5FOBRZb7XjCCiKedruLwCgwJ3f.png) ![20ts|410x500](upload://fe2VtpRvSHUq7BH2GBBwN7kaek1.png)
4

The initial condition is to power(-1). I have also attached the results again. Sorry for the wrong email.

1
1.png554×674 4.29 KB

5

I Googled the problem and found this paper:

https://aip.scitation.org/doi/10.1063/1.434833

The following short example based on that paper seems to have the desired qualitative behavior:

from dolfin import *

mesh = UnitSquareMesh(100,100)
V = FunctionSpace(mesh,"CG",1)
n = Function(V)
n_old = Function(V)
n_old.rename("n","n")
dn = TestFunction(V)

C = Constant(1.0)
kappa = Constant(0.05**2)
n2 = Constant(0.25)
n3 = Constant(1.25)
df0_dn = 4.0*C*n*(n-n2)*(n-n3)

dF = (df0_dn*dn + 2.0*kappa*inner(grad(n),grad(dn)))*dx

Dt = Constant(1.0e-1)
dn_dt = (n-n_old)/Dt

res = inner(dn_dt,dn)*dx + dF

x = SpatialCoordinate(mesh)
n_init = n3/(1.0 + exp(n3*sqrt(C/kappa)*(x[0]-0.5)))
n_old.assign(project(n_init,V))
n.assign(n_old)

nFile = File("results/n.pvd")
for i in range(0,100):
    nFile << n_old
    solve(res==0,n)
    n_old.assign(n)

Hope that helps!

1 Like
6

Dear David,
Thanks so much for your help. It actually worked. However, I still did not understand what actually goes wrong in my code. It is still a question for me.

Regards,
Benhour