Newer
Older
# Extended Standard Cahn-Hilliard Example to Binary Flory Huggins
# as discussed on 28/11/2019
# Example can be found at
# https://bitbucket.org/fenics-project/dolfin/src/master/python/demo/
# documented/cahn-hilliard/demo_cahn-hilliard.py.rst#rst-header-id1
# Runs with fenics 2019.01
# The resulting .pvd file can be opened using default settings in ParaView
import matplotlib.pyplot as plt
#import mshr as ms
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
import numpy as np
from scipy.optimize import curve_fit
import time
import random
from dolfin import *
def create_expr(p_list):
string = ''
for p in p_list:
string = (string+'(x[0]-'+str(p[0])+')*(x[0]-'+str(p[0])+')+(x[1]-'+str(p[1])+
')*(x[1]-'+str(p[1])+')+(x[2]-'+str(p[2])+')*(x[2]-'+str(p[2])+
')<=.05*.05 ? .76 :' )
string = string + '.268'
return string
# test = np.loadtxt('/Users/hubatsch/ownCloud/Dropbox_Lars_Christoph/frap_theory/points_clean.txt', delimiter=',')
test = [[0.5, 0.5, 0.5]]
f = Expression((create_expr(test), 1), degree=1)
# Class for interfacing with the Newton solver
class CahnHilliardEquation(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)
# Model parameters
kappa = 0.25*10e-05 # surface parameter
dt = 0.1e-02 # time step
X = 2.2 # Chi Flory Huggins parameter
# time stepping family, e.g.:
# theta=1 -> backward Euler, theta=0.5 -> Crank-Nicolson
theta = 0.5
# Form compiler optionsmai
parameters["form_compiler"]["optimize"] = True
parameters["form_compiler"]["cpp_optimize"] = True
# mesh = Mesh('Meshes/te.xdmf')
mesh = Mesh()
with XDMFFile("/home/hubatsch/frap_theory/Meshes/box_with_sphere.xdmf") as infile:
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
infile.read(mesh)
# domain = ms.Sphere(Point(0, 0, 0), 1.0)
# mesh = ms.generate_mesh(domain, 100)
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
ME = FunctionSpace(mesh, P1*P1)
# Define trial and test functions
du = TrialFunction(ME)
q, v = TestFunctions(ME)
# Define functions
u = Function(ME) # current solution
u0 = Function(ME) # solution from previous converged step
tis = time.time()
# Split mixed functions
dc, dmu = split(du)
c, mu = split(u)
c0, mu0 = split(u0)
# Create intial conditions and interpolate
# u_init = InitialConditions(degree=1)
# u.interpolate(u_init)
# u0.interpolate(u_init)
u.interpolate(f)
u0.interpolate(f)
print(time.time()-tis)
# Compute the chemical potential df/dc
c = variable(c)
# mu_(n+theta)
mu_mid = (1.0-theta)*mu0 + theta*mu
# Weak statement of the equations
L0 = c*q*dx - c0*q*dx + dt*c*(1-c)*dot(grad(mu_mid), grad(q))*dx
L1 = mu*v*dx - (ln(c/(1-c))+X*(1-2*c))*v*dx - kappa*dot(grad(c), grad(v))*dx
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 = CahnHilliardEquation(a, L)
solver = NewtonSolver()
# solver.parameters["linear_solver"] = "lu"
solver.parameters["linear_solver"] = 'gmres'
#solver.parameters["preconditioner"] = 'ilu'
#solver.parameters["convergence_criterion"] = "incremental"
#solver.parameters["relative_tolerance"] = 1e-6
# Output file
# file = File("output.pvd", "compressed")
file_c = XDMFFile('c_long.xdmf')
print(mesh.num_cells())
# Step in time
t = 0.0
T = 3*dt
ti = time.time()
while (t < T):
# file << (u.split()[0], t)
file_c.write(u.split()[0], t)
t += dt
u0.vector()[:] = u.vector()
solver.solve(problem, u.vector())
print(time.time() - ti)
file_c.close()