Hello there, Could someone please help me with my current fenics problem? Actually I’m working on the PDE problem:
v \in (0, 1), \quad u_z: u(x,y), \quad \mathbb{R}^n \supset \tilde{\omega} \supset \omega\quad \int_{\omega} (v^2 + \eta_{\epsilon})|\nabla u_z|^2 \, dx + G_c \int_{\tilde{\omega}} \left(\frac{1-v^2}{4\epsilon} + \epsilon |\nabla v|^2\right) \, dx
This is the strong form of the problem and the first integral is with respect an elastic energy term and the other one is about fracture energy. Applying the Gateaux derivative and simplifying, I got the following system of equations:
\int_{\omega} (v^2 + \eta_{\epsilon}) \nabla u \cdot \nabla \psi_1 \, dx = 0
\int_{\omega} |\nabla u|^2 \cdot v \cdot \psi_2 \, dx + G_c \int_{\omega} \left(v \cdot \psi_2 - \psi_2\right) \cdot \frac{1}{2\epsilon} + 2\epsilon \cdot \nabla \psi_2 \cdot \nabla v \, dx = 0
And I have the restriction that the fracture term cannot decrease, so I have to start with the function v, by definition of the problem as 1 and the next solution, needs to be lesser than 1 (because there is 1-v² in the fracture energy term). The code that I’ve done, is the following:
from __future__ import print_function
from petsc4py import PETSc
from fenics import *
import matplotlib.pyplot as plt
import numpy as np
import time
combined_file = XDMFFile("combined_solution.xdmf")
combined_file.parameters["flush_output"] = True # Ensure data is written immediately
T = 25.0 # final time
num_steps = 500 # number of time steps
dt = T / num_steps # time step size
tol = 1E-14; eta = 1*10E-7; eps = 4*10E-2; delta = 0.5; Alternate_tol = 1*10E-2; max_iter = 40
G_c = 1
mesh = UnitSquareMesh(5,5)
mesh_h = UnitSquareMesh(5,5)
V = FunctionSpace(mesh, "Lagrange", 1)
g_t = Expression('t', degree=1, t=0)
def boundary_D(x, on_boundary):
return on_boundary and (near(x[0], 0, tol) or near(x[0], 5, tol))
bc = DirichletBC(V, g_t, boundary_D)
u = Function(V)
v = Function(V)
u1 = TestFunction(V)
v1 = TestFunction(V)
v_previous = Function(V)
v_previous.interpolate(Constant(1.0)) # Initialize v0 as 1.0
F1 = (v_previous**2 + eta)*u*u1*dx
F2 = ((grad(u)**2)*v*v1 + 2*G_c*(v-1)*v1/(4*eps) - eps*inner(grad(v), grad(v1)))*dx
vtkfile = File('Brittle1.pvd')
for t in range(num_steps):
# Update the projection of u at the beginning of each time step
check_tol = False
i = 0
while i <= max_iter:
maxnorm = 1
solve(F1 == 0, u, bc)
# Solve for v using the previous time step's solution as an upper bound
solve(F2 == 0, v)
u_vec = as_backend_type(u.vector()).vec()
v_vec = as_backend_type(v.vector()).vec()
v_previous_vec = as_backend_type(v_previous.vector()).vec()
diff_vec = v_vec.copy()
diff_vec.axpy(-1.0, v_previous_vec)
maxnorm = diff_vec.norm(PETSc.NormType.INFINITY) # Calculate the maximum norm
print('\n',f"Iteration {i}: MaxNorm(v) = {maxnorm}",'\n')
v_local = v.vector().get_local()
v_previous_local = v_previous.vector().get_local()
condition_1 = np.any(v_local < v_previous_local)
condition_2 = np.all((v_local >= 0) & (v_local <= 1))
if i > 0:
if (condition_1 and condition_2):
v_previous.vector().set_local(np.minimum(v.vector().get_local(), v_previous.vector().get_local())) # Adjust values within the v_previous Function
v_previous.assign(v) # Update v_previous with the values of v
i += 1
if(maxnorm < Alternate_tol):
check_tol = True
break
i = i + 1
# Update for the next time step
t += dt
g_t.t = t
combined_file.write(u, t) # Save solution of u
combined_file.write(v, t) # Save solution of v
vtkfile << (v, t)
combined_file.close()
plot(u)
plot(v)
print(v.vector().get_local(),'\n')
print(u.vector().get_local())
My idea was to use interpolation for setting v as 1, and then, mix Picard iteration method with newton’s one (that is already being done in fenics), and at each time step, iterating until a it gets to a convergence criterion (In my case, the max norm between the v and v _previous) and update v only if it’s not bigger than the actual solution, because fractures are not allowed to help.
But, the problem that I’ve been stuck with for straight days, is that the method doesn’t converge, the error only increases, and I never get the right fracture evolution, and I don’t know why. Could someone help me please with this problem?