Time-dependent Ginzburg-Landau equation

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 = (2
k
Almbda(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

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$.

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)

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

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!

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