An arbitrary function approaches the critical nucleus in the Ginzburg-Landau time dependent equation

critnuc.py

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)