Skip to content

Instantly share code, notes, and snippets.

@noamraph
Created June 4, 2013 12:46
Show Gist options
  • Save noamraph/5705624 to your computer and use it in GitHub Desktop.
Save noamraph/5705624 to your computer and use it in GitHub Desktop.
An arbitrary function approaches the critical nucleus in the Ginzburg-Landau time dependent equation
from __future__ import division
import numpy as np
from scipy.optimize import brentq
from fipy import Grid1D, CellVariable, Viewer, TransientTerm, DiffusionTerm
import matplotlib.pyplot as plt
D=1.
dx = 0.1
nx = 1000
dt = 0.1
mesh = Grid1D(dx=dx, nx=nx)
x = mesh.cellCenters[0] - dx*nx/2
u_var = CellVariable(name='u', mesh=mesh)
u_c_var = CellVariable(name='u_c', mesh=mesh)
viewer = Viewer(vars=[u_var, u_c_var], datamin=-1, datamax=1)
viewer.limits['datamax'] = 1.5
viewer.limits['datamin'] = -0.5
eq = TransientTerm() == DiffusionTerm(coeff=D) - u_var*(u_var-1/4)*(u_var-1)
u_star = 1/6*(5-np.sqrt(7))
get_x_c = lambda u: np.sqrt(D/2)*np.sqrt(2)*(np.log(7)+2*np.log(u)-2*np.log(3-5*u+np.sqrt(9+6*u*(-5+3*u))))
get_u_c = lambda x: brentq(lambda u: get_x_c(u)+abs(x), 0, u_star) if x != 0 else u_star
u_c = np.array(map(get_u_c, x))
u_c_var.setValue(u_c)
gaussian = np.exp(-x**2)
assymetric = np.choose(x > 0, [np.exp(-x**2), np.exp(-(x/10)**2)])
get_V = lambda u: 1/8*u**2 - 5/12*u**3 + 1/4*u**4
diff = lambda u: (np.roll(u, -1)-np.roll(u, 1))/(2*dt)
get_F = lambda u: np.sum(get_V(u) + D/2*diff(u)**2)
def closest_uc(ar):
"""Return the closest critical nucleus to an array"""
return np.roll(u_c, np.convolve(ar, u_c).argmax()+1)
def solve_initial(initial, decide_fast=True):
history = [initial.copy()]
u_var.setValue(initial)
while u_var.value[450:550].max()-u_var.value[450:550].min() > 0.01:
eq.solve(var=u_var, dt=dt)
history.append(u_var.value.copy())
u_c_var.setValue(closest_uc(u_var.value))
viewer.plot()
if decide_fast:
if get_F(u_var.value) < 0.37:
break
return np.array(history)
def binsearch_step(bottom, top, template):
if top is not None:
A = (top+bottom)/2
else:
A = bottom*2 if bottom > 0 else 1
print A
history = solve_initial(template*A)
#is_high = history[-1][500] > 0.5
is_high = max(history[-1]) > u_star
if is_high:
top = A
else:
bottom = A
return bottom, top, A, history
def binsearch(bottom=0, top=None, n=10, template=assymetric):
As = []
histories = []
for _i in range(n):
bottom, top, A, history = binsearch_step(bottom, top, template)
As.append(A)
histories.append(history)
return bottom, top, As, histories
def graphs(history):
T = np.arange(len(history)) * dt
slowdown = np.abs(history[1:] - history[:-1]).max(axis=1) * dt
slowdown_fig = plt.figure()
plt.gca().set_yscale('log')
plt.plot(T[:len(slowdown)], slowdown)
plt.xlabel('$t$')
plt.ylabel(r'$\max_x\|\frac{\partial u}{\partial t}\|$',
rotation='horizontal', fontsize='large')
plt.title('Critial slowdown')
closest_uc_1 = closest_uc(history[np.argmin(slowdown)])
approach_fig = plt.figure()
plt.yticks([])
times = [0]+[2**i for i in range(8)]
for i, t in enumerate(times):
plt.plot(x[300:700], len(times)-1-i + history[t][300:700], 'b-')
plt.plot(x[300:700], len(times)-1-i + closest_uc_1[300:700], 'g-')
plt.text(x[350], len(times)-1-i+0.2, 't=%.1f' % (dt*t))
plt.xlabel('$x$')
plt.ylim(-0.5, len(times))
plt.title('Approaching the critical nucleus')
dist = np.abs(history - closest_uc_1).max(axis=1)
dist_fig = plt.figure()
plt.gca().set_yscale('log')
plt.plot(T, dist)
plt.xlabel('$t$')
plt.ylabel(r'$\max_x\|u-u_c\|$',
rotation='horizontal', fontsize='large')
plt.title('Distance from the critical nucleus')
F = np.array(map(get_F, history))
F_fig = plt.figure()
plt.plot(T, F)
plt.xlabel('$t$')
plt.ylabel('$F$')
plt.title('Decrease in free energy')
return slowdown_fig, approach_fig, dist_fig, F_fig
def main():
bottom, top, _As, _histories = binsearch(n=17)
A = (bottom+top)/2
print A
history = solve_initial(A*assymetric, decide_fast=False)
graphs(history)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment