test_par.py as found from previous work. Paths not functional.

test_par.py

+# 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
+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("/Users/hubatsch/ownCloud/Dropbox_Lars_Christoph/frap_theory/mesh.xdmf") as infile:
+    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()