#code to run the simulation of the system with the Phase Field Model, Diffusion of CU ions and without Diffusion of Additives
import matplotlib.pyplot as plt
from dolfin import *
import numpy as np
import ufl_legacy as ufl
from ufl_legacy import nabla_div
from openpyxl import Workbook
from mshr import *
from tqdm import tqdm
import os
#random initial condition generator
def generate_random_initial(n):
return 0.5np.ones(n) + 0.02np.random.random_sample(n) - 0.01
def G(xi: float) → float:
“”“Calculate G function for the phase field model.”“”
return (4xi**3 - 6xi**2 + 2*xi)
def H(xi):
return (30xi**4 -60xi3 +30*xi2)
def a(xb):
return (1)
def b(xb):
return (1-xb)*xb
#for write now i am removing input param reading directly from code.
velocity = 0
length_of_dom = 5
width_of_dom = 1
mesh_size = 50
delt = 1e-3
sim_time = 10
inter_loc = 4
inter_loc_cp = 5
w_fctr = 10
l_sigma = 2e-1
l_eta = 5e-3
alpha = 0.5
kappa = 1e-2
mobility = 1
mobility2 = 1
Wb = -10
Wa = .25
Wb2 = -45
Wa2 = .25
Ds = 1
cinit = 10
xb_t0 = 0.01
xb2_t0 = 0.01
run_number = 177
save_step = 2
interpolation_degree = 2
Allencahn Constants
dt = delt
nFbyRT = Constant(296485/(8.3144958298))
Copper Diffusion Constants
Vm = 7.11 # molar volume of Cu [cc/mol]
cs = 1.0/Vm
c0 = 1.0/1000.0 # 1 M solution as reference
del_eta = -0.1
RT= Constant(1.0)
xb_eq = xb_t0 * (pow(np.e, (-1.0/16.0)(-100)))
xb2_eq = xb2_t0 * (pow(np.e, (-1.0/16.0)(-100)))
print(‘xb_eq’,xb_eq)
Load initial conditions from NPZ file
print(‘######### loading data file #######’)
f = np.load(“initial_conditions.npz”)
print(“Available arrays:”, f.files)
print(“dp shape:”, f[‘dp’].shape)
print(f[‘dp’])
print(“data shape:”, f[‘data’].shape)
print(f[‘data’])
print(‘######## data store information ########’)
npz file’data’ contains the initial conditions for the phase-field variables
initial_conditions = f[‘data’] # Shape should be (ng, nx, ny)
Get domain parameters from NPZ file
lx = f[‘dp’][0]
ly = f[‘dp’][1]
nx = int(f[‘dp’][2])
ny = int(f[‘dp’][3])
ng = int(f[‘dp’][8]) # Number of grains
print for dbug
print(f[‘data’].shape)
print(f[‘xg’].shape)
print(f[‘yg’].shape)
print(‘lx=’,lx)
print(‘ly=’,ly)
print(‘nx=’,nx)
print(‘ny=’,ny)
print(‘ng=’,ng)
Check the shape of the initial conditions
print(“Initial conditions shape:”, initial_conditions.shape)
Mesh and Function Space Setup:
Create mesh
mesh = RectangleMesh(Point(0.0, 0.0), Point(lx, ly), nx, ny)
#domain = Rectangle(Point(0.0, 0.0), Point(length_of_dom, width_of_dom))
mesh = RectangleMesh(Point(0.0, 0.0), Point(length_of_dom, width_of_dom), #int(length_of_dommesh_size),int(width_of_dommesh_size))
Function spaces
V = FunctionSpace(mesh, ‘Lagrange’, 2) # Using interpolation_degree = 2
W = VectorFunctionSpace(mesh, ‘Lagrange’, 2)
v = Function(W)
Initialize arrays for multiple grains
xi =
delxi =
xi_old =
eta =
Initialize arrays for additives
xb = Function(V)
delxb = TestFunction(V)
xb2 = Function(V)
delxb2 = TestFunction(V)
Initialize additives
xb_old = Function(V)
xb_old = interpolate(Constant(xb_t0), V)
xb2_old = Function(V)
xb2_old = interpolate(Constant(xb2_t0), V)
Get grain IDs from the data
unique_grains = np.unique(f[‘data’])
unique_grains = np.delete(unique_grains, 0) # remove 0 values
grain_indices = {}
Create grain mapping
for idx, grain_id in enumerate(unique_grains):
grain_indices[int(grain_id)] = idx + 1
num_grains = len(grain_indices)
print(“Detected unique_grains:”, unique_grains)
print(“Detected grain mapping:”, grain_indices)
print(“Number of grains:”, num_grains)
#Part 4 - Variable Initialization:
Initialize phase field variables for each grain
for i in range(ng):
xi.append(Function(V))
delxi.append(TestFunction(V))
xi_old.append(Function(V))
eta.append(Function(V))
# Get coordinates of mesh vertices
dof_coordinates = V.tabulate_dof_coordinates()
values = np.zeros(V.dim())
# Initialize values for each grain
for j, coord in enumerate(dof_coordinates):
x_idx = int(coord[0] / lx * (nx - 1))
y_idx = int(coord[1] / ly * (ny - 1))
x_idx = min(max(x_idx, 0), nx-1)
y_idx = min(max(y_idx, 0), ny-1)
grain_index = f['data'][y_idx, x_idx]
if grain_index in grain_indices:
if grain_indices[grain_index] == (i + 1):
values[j] = 1
else:
values[j] = 0
else:
values[j] = 0
xi_old[i].vector()[:] = values
eta[i] = del_eta * xi_old[i]
# Print some debug information
print(f"Grain {i}:")
print(f"Min value: {values.min()}")
print(f"Max value: {values.max()}")
print(f"Number of DOFs: {V.dim()}")
# Verify the initialization
print(f"Grain {i} norm: {xi_old[i].vector().norm('l2')}")
boundary_marker = MeshFunction(‘size_t’, mesh, mesh.topology().dim() - 1)
class lbc(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-15
return on_boundary and abs(x[0]) < (tol)
class rbc(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-15
return on_boundary and abs(length_of_dom - x[0]) < (tol)
class tbc1(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-15
return on_boundary and abs(width_of_dom - x[1]) < (tol) and (x[0] - inter_loc) < (tol)
class tbc2(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-15
return on_boundary and abs(width_of_dom - x[1]) < (tol) and (x[0] - inter_loc) > (tol)
class bbc(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-15
return on_boundary and abs(x[1]) < (tol)
bx0 = rbc()
bx1 = lbc()
bx21 = tbc1()
bx22 = tbc2()
bx3 = bbc()
bx0.mark(boundary_marker,0)
bx1.mark(boundary_marker,1)
bx21.mark(boundary_marker,2)
bx3.mark(boundary_marker,3)
bx22.mark(boundary_marker,4)
ds = Measure(‘ds’, domain = mesh, subdomain_data = boundary_marker)
defining Boundary Conditions for allen cahn equation
boundary_conditions_xi = {
0 : {‘Neumann’: 0.0},
1 : {‘Dirichlet’:(1.0,bx1)},
2 : {‘Dirichlet’:(1.0,bx21)},
3 : {‘Neumann’: 0.0},
4 : {‘Neumann’: 0.0}}
Boundary conditions for each grain
bcsxi =
for grain in range(ng):
grain_bcs =
for i in boundary_conditions_xi:
if ‘Dirichlet’ in boundary_conditions_xi[i]:
expression, sub_domain = boundary_conditions_xi[i][‘Dirichlet’]
bc = DirichletBC(V, expression, sub_domain)
grain_bcs.append(bc)
bcsxi.append(grain_bcs)
integral_N1 =
for grain in range(ng):
grain_integrals =
for i in boundary_conditions_xi:
if ‘Neumann’ in boundary_conditions_xi[i]:
grad_xi = boundary_conditions_xi[i][‘Neumann’]
grain_integrals.append(kappagrad_xidelxi[grain]*ds(i))
integral_N1.append(grain_integrals)
Boundary Conditions For Copper Diffusion
Modify the copper diffusion calculations for multiple grains
Dxi = [Ds * (1-xi_old[i]**4) for i in range(ng)]
m1_xi = [mobility * (1 - xi_old[i]**4) for i in range(ng)]
m2_xi = [mobility2 * (1 - xi_old[i]**4) for i in range(ng)]
Boundary Conditions For Additives
boundary_conditions_xb = {
0 : {‘Dirichlet’: (xb_t0, bx0)},
1 : {‘Neumann’: 0.0},#{‘Dirichlet’: (0.0, bx1)},
2 : {‘Neumann’: 0.0},#{‘Dirichlet’: (xb_t0, bx21)},
3 : {‘Neumann’: 0.0},#{‘Dirichlet’: (xb_t0, bx3)},
4 : {‘Neumann’: 0.0}}#{‘Dirichlet’: (xb_t0, bx22)}}
bcsxb =
for i in boundary_conditions_xb:
if ‘Dirichlet’ in boundary_conditions_xb[i]:
expression, sub_domain = boundary_conditions_xb[i][‘Dirichlet’]
bc = DirichletBC(V, expression, sub_domain)
bcsxb.append(bc)
integral_N2 =
integral_N3 =
for i in boundary_conditions_xb:
if ‘Neumann’ in boundary_conditions_xb[i]:
grad_xb = boundary_conditions_xb[i][‘Neumann’]
# Create separate integrals for each grain
for j in range(ng):
integral_N2.append(a(xb_old)m1_xi[j]grad_xbdelxbds(i))
integral_N3.append(b(xb_old)*m1_xi[j]G(xi[j])grad_xbdelxbds(i))
Boundary Conditions For Additives2
boundary_conditions_xb2 = {
0 : {‘Dirichlet’: (xb2_t0, bx0)},
1 : {‘Neumann’: 0.0},#{‘Dirichlet’: (0.0, bx1)},
2 : {‘Neumann’: 0.0}, #{‘Dirichlet’: (xb2_t0, bx21)},
3 : {‘Neumann’: 0.0}, #{‘Dirichlet’: (xb2_t0, bx3)},
4 : {‘Neumann’: 0.0}} #{‘Dirichlet’: (xb2_t0, bx22)}}
bcsxb2 =
for i in boundary_conditions_xb2:
if ‘Dirichlet’ in boundary_conditions_xb2[i]:
expression, sub_domain = boundary_conditions_xb2[i][‘Dirichlet’]
bc = DirichletBC(V, expression, sub_domain)
bcsxb2.append(bc)
integral_N4 =
integral_N5 =
for i in boundary_conditions_xb2:
if ‘Neumann’ in boundary_conditions_xb2[i]:
grad_xb2 = boundary_conditions_xb2[i][‘Neumann’]
# Create separate integrals for each grain
for j in range(ng):
integral_N4.append(a(xb2_old)m2_xi[j]grad_xb2delxb2ds(i))
integral_N5.append(b(xb2_old)*m2_xi[j]G(xi[j])grad_xb2delxb2ds(i))
#Weak Form of Allen Cahn equation, Additives, Copper Diffusion
Modify the Fxi definition to include proper grain interactions
Fxi =
for i in range(ng):
other_grains = sum(xi_old[j] for j in range(ng) if j != i)
F = (
xi[i] * delxi[i] * dx - xi_old[i] * delxi[i] * dx + dt * l_sigma * (Wb2*xb2_old +
Wb*xb_old + Wa*(1 - xb_old - xb2_old)) * G(xi[i]) * delxi[i] * dx
+ dt * l_sigma * inner(nabla_grad(delxi[i]), kappa*nabla_grad(xi[i])) * dx
- dt * l_sigma * sum(integral_N1[i]) + dt * l_eta * (10 + (100*(xb - xb2))) *
H(xi_old[i]) * ((np.e**((1-alpha)*nFbyRT*(eta[i]))) - (cinit)*np.e**(-
alpha*nFbyRT*(eta[i]))) * delxi[i] * dx + dt * Constant(10.0) * xi[i] *
other_grains * delxi[i] * dx # Increased interaction strength
)
Fxi.append(F)
Fxb = (
xb * delxb * dx - xb_old * delxb * dx - dt * RT * sum(integral_N2) + dt *
sum(m1_xi[i] * RT * inner(nabla_grad(delxb), a(xb_old)*nabla_grad(xb)) for i in
range(ng)) * dx - dt * (Wb - Wa) * sum(integral_N3) + dt * sum(m1_xi[i] * (Wb - Wa) *
inner(nabla_grad(delxb), b(xb_old)*G(xi[i])*nabla_grad(xi[i])) for i in range(ng)) *
dx)
Fxb2 = (
xb * delxb2 * dx - xb2_old * delxb2 * dx - dt * RT * sum(integral_N4) + dt *
sum(m2_xi[i] * RT * inner(nabla_grad(delxb2), a(xb2_old)*nabla_grad(xb2)) for i in
range(ng)) * dx - dt * (Wb2 - Wa2) * sum(integral_N5) + dt * sum(m2_xi[i] * (Wb2 -
Wa2) *inner(nabla_grad(delxb2), b(xb2_old)*G(xi[i])*nabla_grad(xi[i])) for i in
range(ng))
*dx)
#Fxb = (
xb * delxb * dx
- xb_old * delxb * dx
+ dt * mobility * RT * inner(nabla_grad(delxb), a(xb_old)*nabla_grad(xb)) * dx
+ dt * mobility * (Wb - Wa) * inner(nabla_grad(delxb), b(xb_old)*sum(G(xi[i]) for i in range(ng))*nabla_grad(sum(xi[i] for i in range(ng)))) * dx
#)
#Fxb2 = (
xb2 * delxb2 * dx
- xb2_old * delxb2 * dx
+ dt * mobility2 * RT * inner(nabla_grad(delxb2), a(xb2_old)*nabla_grad(xb2)) * dx
+ dt * mobility2 * (Wb2 - Wa2) * inner(nabla_grad(delxb2), b(xb2_old)*sum(G(xi[i]) for i in range(ng))*nabla_grad(sum(xi[i] for i in range(ng)))) * dx
#)
solver_parameters = {
“newton_solver”: {
“relative_tolerance”: 1e-5,
“absolute_tolerance”: 1e-6,
“maximum_iterations”: 100,
“relaxation_parameter”: 0.5,
“linear_solver”: “mumps”,
“krylov_solver”: {
“absolute_tolerance”: 1e-7,
“relative_tolerance”: 1e-6,
“maximum_iterations”: 1000
}
}
}
Create output directory
output_dir = f’trial1/study_number{run_number}’
os.makedirs(output_dir, exist_ok=True)
Output files setup
xi_file = File(f’{output_dir}/order_parameter.pvd’)
acc_file = File(f’{output_dir}/accelerator.pvd’)
supp_file = File(f’{output_dir}/suppressor.pvd’)
Time evolution loop
set_log_level(1000)
iter = 0
t = 0.0
for t in np.arange(0, sim_time, delt):
# Solve for each grain
try:
for i in range(ng):
solve(Fxi[i] == 0, xi[i], bcsxi[i], solver_parameters=solver_parameters)
solve(Fxi[i] == 0, xi[i], bcsxi[i])
# Solve for additives
solve(Fxb == 0, xb, bcsxb)
solve(Fxb2 == 0, xb2, bcsxb2)
solve(Fxb == 0, xb, bcsxb, solver_parameters=solver_parameters)
solve(Fxb2 == 0, xb2, bcsxb2, solver_parameters=solver_parameters)
except RuntimeError as e:
print(f"Solver failed at t={t}: {e}")
break
if iter % save_step == 0:
# Calculate combined field for visualization
xi_tot = Function(V, name='Combined_Grains')
xi_tot.assign(sum((2*i + 1) * xi_old[i] for i in range(ng)))
xi_file << (xi_tot, t)
print(f"Time step: {t:.3f}/{sim_time}")
to save all xi values in one file
for i in range(ng):
xi_func = Function(V, name=f'xi{i+1}')
xi_func.assign(xi[i])
xi_file << (xi_func, t)
acc_file << (xb_old, t)
supp_file << (xb2_old, t)
print('writing ', t)
# Update old values
for i in range(ng):
xi_old[i].assign(xi[i])
xb_old.assign(xb)
xb2_old.assign(xb2)
iter += 1
Final output
xi_tot = Function(V, name=‘Combined_Grains’)
xi_tot.assign(sum((2*i + 1) * xi_old[i] for i in range(ng)))
xi_file << (xi_tot, t)
acc_file << (xb_old, t)
supp_file << (xb2_old, t)