Hi all,
Any guidance for this error message? I am a pretty novice coder, and am a student researcher working with a professor using the FEM on the Phase Field Crystal Equation. I am more than happy to send my code if that would help.
One would definitely need the code to give any further guidance.
Please modify your post and format code with
```python
# add code here
```
#Create mesh
msh = mesh.create_rectangle(
MPI.COMM_WORLD,
points=[(0.0, 0.0), (1.0, 1.0)],
n=(16, 16),
cell_type=CellType.triangle,
ghost_mode=GhostMode.shared_facet,
)
tdim = msh.topology.dim
msh.topology.create_connectivity(tdim - 1, tdim)
#facets = msh.exterior_facet_indices(msh.topology)
P2 = element(“Lagrange”, msh.basix_cell(), 2)
ME = functionspace(msh, mixed_element([P2, P2]))
xdmf = io.XDMFFile(msh.comm, “diffusion.xdmf”, “w”)
xdmf.write_mesh(msh)
#Define trial and test functions
dc, dmu = ufl.TrialFunctions(ME)
q, v = ufl.TestFunctions(ME)
#Define functions
u = Function(ME) # current solution u=[c;mu]
du = Function(ME)
u0 = Function(ME) # solution from previous converged step
u_newt = Function(ME) # solution from previous Newton step
#Split mixed functions
c, mu = ufl.split(u)
c0, mu0 = ufl.split(u0)
c_newt, mu_newt = ufl.split(u_newt)
#Define exact initial conditions
class exact_solution0:
def __ call __(self, x):
return
class exact_solution1:
def __ call __(self, x):
return these equations are very complex, excluded for brevity
u_exact0 = exact_solution0()
u_exact1 = exact_solution1()
#Set initial condition
u.sub(0).interpolate(u_exact0)
u0.sub(0).interpolate(u_exact0)
u_newt.sub(0).interpolate(u_exact0)
u.sub(1).interpolate(u_exact1)
u0.sub(1).interpolate(u_exact1)
u_newt.sub(1).interpolate(u_exact1)
u.x.scatter_forward()
u0.x.scatter_forward()
u_newt.x.scatter_forward()
xdmf.write_function(u_newt, 0)
h = CellDiameter(mesh)
h_avg = ((h(‘+’) + h(‘-’))/2.0)
n = FacetNormal(mesh)
alpha = 10.0
#Define variational forms
aIP_LHS = (
inner(div(grad(dc)), div(grad(v)))*dx
- inner(avg(div(grad(dc))), jump(grad(v), n))*dS
- inner(jump(grad(dc), n), avg(div(grad(v))))*dS
+ (alpha/h_avg)*inner(jump(grad(dc), n), jump(grad(v), n))*dS
- inner(div(grad(dc)), dot(grad(v), n))*ds
- inner(dot(grad(dc), n), div(grad(v)))*ds
+ (alpha/h)*inner(dot(grad(dc), n), dot(grad(v), n))*ds
)
aIP_RHS = (
inner(div(grad(c_newt)), div(grad(v)))*dx
- inner(avg(div(grad(c_newt))), jump(grad(v), n))*dS
- inner(jump(grad(c_newt), n), avg(div(grad(v))))*dS
+ (alpha/h_avg)*inner(jump(grad(c_newt), n), jump(grad(v), n))*dS
- inner(div(grad(c_newt)), dot(grad(v), n))*ds
- inner(dot(grad(c_newt), n), div(grad(v)))*ds
+ (alpha/h)*inner(dot(grad(c_newt), n), dot(grad(v), n))*ds
)
#Weak statement of the equations for LHS and RHS
L0 = dc*q*dx + dt*dot(grad(dmu), grad(q))*dx
L1 = dmu*v*dx - 3*c_newt**2dc*v*dx - (1-eps)*dc*v*dx - aIP_LHS
L = L0 + L1
f0 = c_newt*q*dx - c0*q*dx + dt*dot(grad(mu_newt), grad(q))*dx
f1 = mu_newt*v*dx - c_newt**3vdx - (1-eps)*c_newt*v*dx + 2*dot(grad(c0), grad(v))*dx - aIP_RHS
F = f0 + f1
#Solver
problem = LinearProblem(L, F)
res_tol = 1e-6
tot_energy_array = []
newton loop
t = 0 # Initial Time
T = 100 # Final Time
Nsteps = 100 # Number of Time Steps
dt = T / Nsteps # Time Step Size
i = 0
while t < T:
i += 1
t += dt
res_err = 1
k = 0
while res_err > res_tol:
k += 1
du = problem.solve()
u_newt.x.array[:] -= du.x.array[:]
res_err = la.norm(u_newt.x.array - u.x.array) / la.norm(u.x.array)
print(f"Residual error: {res_err}, Iteration: {k}")
u.x.array[:] = u_newt.x.array[:]
u0.x.array[:] = u.x.array[:]
tot_energy_array.append(t)
print(f"Step: {i}, Time: {t}")
xdmf.write_function(u_newt, t)
xdmf.close()