"""This is the solver in support of the paper "Phase field model for coupled displacive and diffusive microstructural processes under thermal loading" The code was run with the development versions of UFC, UFC, FFC and DOLFIN from the FEniCS Project (http://www.fenicsproject.org) in October 2010. """ __author__ = "Garth N. Wells (gnw20@cam.ac.uk) and Mirko Maraldi " __date__ = "2010-10-09" __copyright__ = "Copyright (C) 2010 Garth N. Wells and Mirko Maraldi" __license__ = "GNU LGPL Version 2.1" import os, random from dolfin import * # Some FFC parameters parameters["form_compiler"]["log_level"] = INFO parameters["form_compiler"]["cpp_optimize"] = True parameters["form_compiler"]["cpp_optimize_flags"] = "-O3" parameters["form_compiler"]["optimize"] = True # Use PETSc as linear algebra backend parameters["linear_algebra_backend"] = "PETSc" def update(U, U0, v0, a0, beta, gamma, dt): """Update fields at the end of each time step.""" # Get sub-functions (deep copy) u, u0 = U.split(True)[0], U0.split(True)[0] # Get vectors (references) u_vec, u0_vec = u.vector(), u0.vector() v0_vec, a0_vec = v0.vector(), a0.vector() # Update acceleration and velocity a_vec = (1.0/(2.0*beta))*( (u_vec - u0_vec - v0_vec*dt)/(0.5*dt*dt) - (1.0-2.0*beta)*a0_vec ) v_vec = dt*((1.0-gamma)*a0_vec + gamma*a_vec) + v0_vec # Update (U0 <- U) v0.vector()[:], a0.vector()[:] = v_vec, a_vec U0.vector()[:] = U.vector() class InitialConditions(Expression): """Initial conditions for each field.""" def __init__(self, cell_size, temperature): random.seed(2 + MPI.process_number()) self.cell_size = cell_size self.temperature = temperature def eval(self, values, x): values[0]= random.uniform(-1.0e-4, 1.0e-4)*self.cell_size # x displacement values[1]= random.uniform(-1.0e-4, 1.0e-4)*self.cell_size # y displacement values[2]= random.uniform(-1.0e-4, 1.0e-4) # chemical concentration values[3]= 0.0 # chemical potential values[4]= self.temperature # temperature def value_shape(self): return (5,) class ModelEquation(NonlinearProblem): def __init__(self, a, L, bcs, newton_solver): NonlinearProblem.__init__(self) self.L = L self.a = a self.bcs = bcs self.newton_solver = newton_solver self.reset_sparsity = True self.ffc_opt = {"quadrature_degree": 2, "representation": "quadrature"} def F(self, b, x): assemble(self.L, tensor=b, form_compiler_parameters=self.ffc_opt) for bc in bcs: bc.apply(b, x) def J(self, A, x): if self.newton_solver.iteration() > 0: self.newton_solver.linear_solver().parameters["preconditioner"]["reuse"] = True assemble(self.a, tensor=A, reset_sparsity=self.reset_sparsity, form_compiler_parameters=self.ffc_opt) for bc in bcs: bc.apply(A) self.reset_sparsity = False #------------------------------------------------------------------------------ # Strain-related order parameter functions # (eps_xx + eps_yy) def e1_f(r): return r[0].dx(0) + r[1].dx(1) # (eps_xx - eps_yy) def e2_f(r): return r[0].dx(0) - r[1].dx(1) # (1/2)(dv/dx + du/dy) def e3_f(r): return 0.5*(r[0].dx(1) + r[1].dx(0)) #------------------------------------------------------------------------------ # Free-energy functions # Chemical free energy def f_diff(c, T, Dc, d): return d['A']*(c**4)/4.0 + d['B']*(T-d['Tp'])*(c**2)/(2.0*d['Tp']) \ + (d['lambda_c']/2.0)*dot(Dc, Dc) # Displacive/elastic free energy def f_disp(e1, e2, e3, De2, c, T, d): return (d['D']/6.0)*(e2**6) - (d['E']/4.0)*(e2**4) \ + (d['F']/2.0)*((T-d['TM'])/d['TM'])*(e2**2) \ + (d['G']/2.0)*e1*(e1- (d['alpha']*(T-d['T_ref']) + d['x1c']*c + d['x12']*e2**2)) \ + (d['H']/2.0)*e3**2 \ + (d['lambda_e']/2.0)*dot(De2, De2) # Thermal free energy def f_therm(T, d): return - (d['cT']/(2.0*d['T_ref']))*(T-d['T_ref'])**2 # Coupling free energy def f_cpl(c, e2, d): return d['x2c']*(c**2)*(e2**2) #------------------------------------------------------------------------------ # Create linear solver and Newton solver # Create PETSc ParaSails preconditioner #pc = PETScPreconditioner("hypre_parasails") #pc.parameters["hypre"]["parasails"]["threshold"] = 0.15 #pc.parameters["hypre"]["parasails"]["levels"] = 0 pc = PETScPreconditioner("ilu") pc.parameters["ilu"]["fill_level"]= 1 # Create Krylov linear solvers and set parameters krylov_solver = PETScKrylovSolver("bicgstab", pc) krylov_solver.parameters["relative_tolerance"] = 1.0e-6 krylov_solver.parameters["error_on_nonconvergence"] = False krylov_solver.parameters["maximum_iterations"] = 250; # Create Newton solver newton_solver = NewtonSolver(krylov_solver, PETScFactory.instance()) newton_solver.parameters["relative_tolerance"] = 1.0e-6 newton_solver.parameters["absolute_tolerance"] = 1.0e-14 newton_solver.parameters["maximum_iterations"] = 10 newton_solver.parameters["error_on_nonconvergence"] = False #------------------------------------------------------------------------------ # Mesh and problem data # Parameters which define each problem case6 = dict(tag="case_6", dt0=1.0e-2, num_steps=10000, \ T_ref=1.2, T_ext=0.3, delta = 1.0e-4, num_cells=64) case7 = dict(tag="case_7", dt0=2.5e-4, num_steps=10000, T_ref=1.2, T_ext=0.3, delta = 5.0e-2, num_cells=128) case8 = dict(tag="case_8", dt0=5.0e-4, num_steps=10000, \ T_ref=1.2, T_ext=0.3, delta = 5.0e-3, num_cells=128) case9 = dict(tag="case_9", dt0=2.0e-3, num_steps=10000, \ T_ref=1.2, T_ext=0.3, delta = 1.0e-3, num_cells=128) case10 = dict(tag="case_10", dt0=5.0e-4, num_steps=10000, \ T_ref=1.2, T_ext=0.3, delta = 2.0e-3, num_cells=64) case11 = dict(tag="case_11", dt0=2.0e-3, num_steps=1100, \ T_ref=1.2, T_ext=0.3, delta = 1.5e-3, num_cells=128) case12 = dict(tag="case_12", dt0=2.0e-2, num_steps=10000, \ T_ref=1.2, T_ext=0.3, delta = 5.0e-4, num_cells=64) case13 = dict(tag="case_13", dt0=5.0e-4, num_steps=10000, \ T_ref=1.2, T_ext=0.45, delta=5.0e-2, num_cells=128) case14 = dict(tag="case_14", dt0=5.0e-2, num_steps=30000, \ T_ref=1.2, T_ext=0.4, delta=5.0e-2, num_cells=128) case15 = dict(tag="case_15", dt0=2.0e-3, num_steps=30000, \ T_ref=1.2, T_ext=0.1, delta=5.0e-2, num_cells=128) case_data = case6 linear_solver = krylov_solver # Restart data restart = False restart_number = 0 dt = Constant(case_data["dt0"]) t = 0.0 step_number = 0 # Free-energy parameters f_data = dict(A=7.3e-3, B=6.6e-3,\ D=0.3, E=0.0031, F=0.0025, \ G=0.05, H=0.05, \ lambda_e=1.0e-8, lambda_c=1.0e-7, \ TM = 0.495, Tp=1.0, T_ref=case_data["T_ref"], cT=2.0e-3, alpha=1.0e-2, \ x12=1.0e-1, x1c=5.0e-3, x2c=8.0e-2) print "Free energy data:" print f_data print "Case data:" print case_data # Thermal parameters kT = 2.0e-2 # thermal conductivity delta = case_data['delta'] # heat boundary flux coefficient T_ext = case_data['T_ext'] # external temperature # Mechanical parameters rho = 5.0e-7 # solid mass density eta = 5.0e-8 # viscous stress parameter # Mass diffusion parameters Mob0 = 1.0e4 # solute mobility Q = 5.0 # mobility parameter # Time stepping parameters rho_infty = 0.7 alpha_m = (2.0*rho_infty -1.0)/(rho_infty+ 1.0) alpha_f = rho_infty/(rho_infty+ 1.0) beta = 0.25*(1.0 - alpha_m + alpha_f)**2 gamma = 0.5 - alpha_m + alpha_f # Mesh and mesh cell related quantities num_cells = case_data["num_cells"] mesh = UnitSquare(num_cells, num_cells) n = FacetNormal(mesh) h = triangle.circumradius*2.0 h_avg = (h('+') + h('-'))/2.0 # Define Dirichlet boundary def boundary(x, on_boundary): return x[0] < (1.0 - 1.0e-6) and x[1] < (1.0 - 1.0e-6) and on_boundary #------------------------------------------------------------------------------ # Define function spaces P2v = VectorFunctionSpace(mesh, "Lagrange", 2) P1 = FunctionSpace(mesh, "CG", 1) ME = MixedFunctionSpace([P2v, P1, P1, P1]) # Define test and trial functions (w, q, r, z) = TestFunctions(ME) dU = TrialFunction(ME) # Current solution and solution from previous converged step U, U0 = Function(ME), Function(ME) # Acceleration and velocity from previous converged step a0, v0 = Function(P2v), Function(P2v) # Split mixed functions u, c, mu, T = split(U) u0, c0, mu0, T0 = split(U0) #------------------------------------------------------------------------------ # Balance of linear momentum (using the generalised-alpha method) # Acceleration and velocity (current) a = (u - u0 - dt*v0 - (0.5-beta)*dt*dt*a0)/(beta*dt*dt) v = v0 + dt*((1.0-gamma)*a0 + gamma*a) # Fields at 'alpha' midpoint (generalised-alpha method) v_alpha, u_alpha, c_alpha, mu_alpha, T_alpha \ = [(1.0-alpha_f)*y + alpha_f*y0 for y, y0 in zip([v, u, c, mu, T], [v0, u0, c0, mu0, T0])] a_alpha = (1.0-alpha_m)*a + alpha_m*a0 # Variables at 'alpha' mid-point e1_alpha, e2_alpha, e3_alpha = variable(e1_f(u_alpha)), variable(e2_f(u_alpha)), variable(e3_f(u_alpha)) De2_alpha = variable(grad(e2_f(u_alpha))) c_alpha, Dc_alpha = variable(c_alpha), variable(grad(c_alpha)) T_alpha = variable(T_alpha) # Free-energy at the 'alpha' mid-point f_alpha = f_disp(e1_alpha, e2_alpha, e3_alpha, De2_alpha, c_alpha, T_alpha, f_data) \ + f_diff(c_alpha, T_alpha, Dc_alpha, f_data) \ + f_therm(T_alpha, f_data) \ + f_cpl(c_alpha, e2_alpha, f_data) # Gradient stress (more precisely, div(Sigma)) Sigma_alpha = diff(f_alpha, De2_alpha) # Linear form L_u = rho*dot(w, a_alpha)*dx L_u += e2_f(w)*eta*e2_f(v_alpha)*dx L_u += e1_f(w)*diff(f_alpha, e1_alpha)*dx \ + e2_f(w)*diff(f_alpha, e2_alpha)*dx \ + 2.0*e3_f(w)*diff(f_alpha, e3_alpha)*dx L_u += dot(grad(e2_f(w)), Sigma_alpha)*dx \ - dot(jump(e2_f(w), n), avg(Sigma_alpha))*dS \ - dot(avg(f_data['lambda_e']*grad(e2_f(w))), jump(e2_f(u_alpha), n))*dS \ + 8.0*(f_data['lambda_e']/h_avg)*dot(jump(e2_f(w)), jump(e2_f(u_alpha)))*dS #------------------------------------------------------------------------------ # Mass diffusion (using the Crank-Nicolson method) # Mid-point values mu_mid, T_mid = variable(0.5*(mu + mu0)), variable(0.5*(T + T0)) # Variables at end point e1, e2, e3 = variable(e1_f(u)), variable(e2_f(u)), variable(e3_f(u)) De2 = variable(grad(e2_f(u))) c, Dc = variable(c), variable(grad(c)) # Temperature-dependent mobility Mob_mid = Mob0*exp(-Q/T_mid) # Free-energy at end-point (t_(n+1)) f = f_disp(e1, e2, e3, De2, c, T, f_data) \ + f_diff(c, T, Dc, f_data) \ + f_therm(T, f_data) \ + f_cpl(c, e2, f_data) # Weak form L_c = q*((c-c0)/dt)*dx + Mob_mid*dot(grad(q), grad(mu_mid))*dx L_mu = r*mu*dx - r*diff(f, c)*dx - dot(grad(r), diff(f, Dc))*dx #------------------------------------------------------------------------------- # Heat equation (using the Crank-Nicolson method) # Fields at a mid-point (Crank-Nicolson method) u_mid, c_mid, mu_mid, T_mid = 0.5*(u + u0), 0.5*(c + c0), 0.5*(mu + mu0), 0.5*(T + T0) # Variables at mid-point e1_mid, e2_mid, e3_mid = variable(e1_f(u_mid)), variable(e2_f(u_mid)), variable(e3_f(u_mid)) De2_mid = variable(grad(e2_f(u_mid))) c_mid, Dc_mid = variable(c_mid), variable(grad(c_mid)) T_mid = variable(T_mid) # Free-energy at mid-point f_mid = f_disp(e1_mid, e2_mid, e3_mid, De2_mid, c_mid, T_mid, f_data) \ + f_diff(c_mid, T_mid, Dc_mid, f_data) \ + f_therm(T_mid, f_data) \ + f_cpl(c_mid, e2_mid, f_data) # Heat flux (Fourier law) q_mid = -kT*grad(T_mid) sigma1_mid, sigma2_mid, sigma3_mid = diff(f_mid, e1_mid), diff(f_mid, e2_mid), diff(f_mid, e3_mid) vsigma2_mid = eta*e2_f((u-u0)/dt) # (Dsigma/DT)(de/dt) D_sigma_DTde_mid = diff(sigma1_mid, T_mid)*e1_f((u-u0)/dt) \ + diff(sigma2_mid, T_mid)*e2_f((u-u0)/dt) \ + 2.0*diff(sigma3_mid, T_mid)*e3_f((u-u0)/dt) # (Dmu/DT)(dc/dt) D_mu_DTdc_mid = diff(diff(f_mid, c_mid), T_mid)*(c-c0)/dt # ds/dT at mid-point (s = -df/dT) ds_dT_mid = -diff(diff(f_mid, T_mid), T_mid) # Weak form L_T = z*T_mid*ds_dT_mid*((T-T0)/dt)*dx \ - z*T_mid*D_sigma_DTde_mid*dx \ - z*T_mid*D_mu_DTdc_mid*dx \ - z*Mob_mid*dot(grad(mu_mid), grad(mu_mid))*dx \ - z*vsigma2_mid*e2_f((u-u0)/dt)*dx \ - dot(grad(z), q_mid)*dx \ + z*delta*(T_mid - T_ext)*ds #------------------------------------------------------------------------------ # Total linear form and Jacobian L = L_u + L_c + L_mu + L_T a = derivative(L, U, dU) #------------------------------------------------------------------------------ # Output files for postprocessing prefix = "./results/" + case_data["tag"] + "_" + str(restart_number) file_c = File(prefix + "/c.pvd", "compressed") file_e2 = File(prefix + "/e2.pvd", "compressed") file_e1 = File(prefix + "/eps_v.pvd", "compressed") file_u = File(prefix + "/u.pvd", "compressed") file_T = File(prefix+ "/temperature.pvd", "compressed") #------------------------------------------------------------------------------ # Start from t = 0 or restart if not restart: # Create intial conditions U_init = InitialConditions(1.0/case_data["num_cells"], case_data["T_ref"]) # Interpolate initial conditions in a finite element space U.interpolate(U_init) U0.interpolate(U_init) # Update kinematic fields. This will cause problems if the initial # conditions are not compatible with the Dirichlet BCs. update(U, U0, v0, a0, float(beta), float(gamma), float(dt)) else: prefix = "./restart/" + case_data["tag"] + "/n" + str(step_number) + "_t" + str(t) U.interpolate(Function(ME, prefix + "/U.xml")) U0.interpolate(Function(ME, prefix + "/U0.xml")) v0.interpolate(Function(P2v, prefix + "/v0.xml")) a0.interpolate(Function(P2v, prefix + "/a0.xml")) #------------------------------------------------------------------------------ # Create nonlinear problem and solve # Displacement boundary conditions bcs = [DirichletBC(ME.sub(0), Constant((0.0, 0.0)), boundary)] # Create nonlinear problem problem = ModelEquation(a, L, bcs, newton_solver) while (step_number < case_data["num_steps"]): t += float(dt) step_number += 1 print "Time, dt, step number:", t, float(dt), step_number tic() num_iterations, converged = newton_solver.solve(problem, U.vector()) print "CPU time for step:", toc() # Reset preconditioner periodically if (step_number % 20) == 0: newton_solver.linear_solver().parameters["preconditioner"]["reuse"] = False newton_solver.linear_solver().parameters["preconditioner"]["same_nonzero_pattern"] = True else: newton_solver.linear_solver().parameters["preconditioner"]["reuse"] = True tic() update(U, U0, v0, a0, float(beta), float(gamma), float(dt)) # Write solution vectors to file for possible restart if step_number % 20 == 0: print "Saving restart data" extension = str(step_number) + "_t" + str(t) + ".xml" prefix = "./restart/" + case_data["tag"] + "/n" + str(step_number) + "_t" + str(t) file_restart_U = File(prefix + "/U.xml") file_restart_U0 = File(prefix + "/U0.xml") file_restart_v0 = File(prefix + "/v0.xml") file_restart_a0 = File(prefix + "/a0.xml") file_restart_U << U.vector() file_restart_U0 << U0.vector() file_restart_a0 << a0.vector() file_restart_v0 << v0.vector() # Write solution to file for visualisation file_e1 << (project(e1_f(U.split()[0]), P1), t) file_e2 << (project(e2_f(U.split()[0]), P1), t) file_u << (U.split()[0], t) file_c << (U.split()[1], t) file_T << (U.split()[3], t) print "End of script.", t, float(dt), final_t