# DOLFIN encountered an error. If you are not able to resolve this issue

**URL:** <https://fenicsproject.discourse.group/t/dolfin-encountered-an-error-if-you-are-not-able-to-resolve-this-issue/16873>\
**Category:** Errors\
**Created:** [January 29, 2025, 12:39pm UTC](https://fenicsproject.discourse.group/t/dolfin-encountered-an-error-if-you-are-not-able-to-resolve-this-issue/16873 "2025-01-29T12:39:09Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![Naveen\_Kumar](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/naveen_kumar/32/8025_2.png) [@Naveen\_Kumar](https://fenicsproject.discourse.group/u/Naveen_Kumar)\
**Post date:** [January 29, 2025, 12:39pm UTC](https://fenicsproject.discourse.group/t/dolfin-encountered-an-error-if-you-are-not-able-to-resolve-this-issue/16873/1 "2025-01-29T12:39:09Z")

</div>

#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.5_np.ones(n) + 0.02_np.random.random\_sample(n) - 0.01

def G(xi: float) → float:  
“”“Calculate G function for the phase field model.”“”  
return (4_xi\*\*3 - 6_xi\*\*2 + 2\*xi)

def H(xi):  
return (30_xi\*\*4 -60_xi **3 +30\*xi** 2)

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(2_96485/(8.3144958_298))

# 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\_dom_mesh\_size),int(width\_of\_dom_mesh\_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(kappa_grad\_xi_delxi[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\_xbdelxb_ds(i))  
integral\_N3.append(b(xb\_old)\*m1\_xi[j]_G(xi[j])grad\_xbdelxb_ds(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\_xb2delxb2_ds(i))  
integral\_N5.append(b(xb2\_old)\*m2\_xi[j]_G(xi[j])grad\_xb2delxb2_ds(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)

---

<div class="post-metadata">

**Author:** ![Stein](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/stein/32/8751_2.png) [@Stein](https://fenicsproject.discourse.group/u/Stein)\
**Post date:** [January 29, 2025, 2:05pm UTC](https://fenicsproject.discourse.group/t/dolfin-encountered-an-error-if-you-are-not-able-to-resolve-this-issue/16873/2 "2025-01-29T14:05:46Z")

</div>

Hi Naveen,

I’m afraid nobody is going to be able to help you with your problem given the current status of your post. Take a look at [Read before posting: How do I get my question answered?](https://fenicsproject.discourse.group/t/read-before-posting-how-do-i-get-my-question-answered/21)

In particular:

- Write a meaningful title
- Phrase a question
- Format your code
- Make your code _minimal_ in reproducibility.
- Write the error message

---

<div class="post-metadata">

**Author:** ![Naveen\_Kumar](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/naveen_kumar/32/8025_2.png) [@Naveen\_Kumar](https://fenicsproject.discourse.group/u/Naveen_Kumar)\
**Post date:** [January 30, 2025, 3:38am UTC](https://fenicsproject.discourse.group/t/dolfin-encountered-an-error-if-you-are-not-able-to-resolve-this-issue/16873/3 "2025-01-30T03:38:36Z")

</div>

thanks for your suggestions… I am new in Fenics. this is my first post. after reading your link i will post again.

---

<div class="post-metadata">

**Author:** ![francesco-ballarin](https://avatars.discourse-cdn.com/v4/letter/f/f05b48/32.png) [@francesco-ballarin](https://fenicsproject.discourse.group/u/francesco-ballarin)\
**Post date:** [January 30, 2025, 7:17am UTC](https://fenicsproject.discourse.group/t/dolfin-encountered-an-error-if-you-are-not-able-to-resolve-this-issue/16873/4 "2025-01-30T07:17:42Z")

</div>

Closed on January 30, 2025, 7:17am UTC