Hello FEniCS community,
I am trying to solve a one-dimensional nonlinear PDE using finite element method. The nonlinear system is solved using the Newton method, but the Newton iterations fail to converge.
My boundary conditions are
y(0,t) = 0,
y_x(1,t) = 0,
y_xx(1,t) = 0.
The initial condition is
y(x,0) = sin(\pi x).
#########
#########
#########
import numpy as np
from fenics import *
import matplotlib.pyplot as plt
nu = Constant(0.1)
def main(T, N, M, degree):
# 1D mesh
mesh = UnitIntervalMesh(N)
# -------------------------------------------------------
# Boundary markers
# -------------------------------------------------------
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
class LeftBoundary(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[0], 0.0)
class RightBoundary(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[0], 1.0)
LeftBoundary().mark(boundaries, 1)
RightBoundary().mark(boundaries, 2)
ds = Measure("ds", domain=mesh, subdomain_data=boundaries)
V = FunctionSpace(mesh, "CG", degree)
bc_left = DirichletBC(V, Constant(0.0), LeftBoundary())
dt = T / M
# Initial condition
u0 = Expression(
"(sin(pi*x[0]))", #sin(pi*x[0])
degree=2
)
uold = interpolate(u0, V)
u = Function(V)
phi = TestFunction(V)
du = TrialFunction(V)
# -------------------------------------------------------
# Coordinates of the degrees of freedom
# -------------------------------------------------------
x = V.tabulate_dof_coordinates().reshape(-1)
idx = np.argsort(x)
x = x[idx]
# -------------------------------------------------------
# Store initial solution (t = 0)
# -------------------------------------------------------
times = [0.0]
U = [uold.vector().get_local()[idx]]
t = 0.0
# -------------------------------------------------------
# Time-stepping
# -------------------------------------------------------
for n in range(M):
t += dt
u.assign(Constant(0.0))
F = (
((u - uold) / dt) * phi * dx
- 0.1 * (
u.dx(0).dx(0)
) * phi.dx(0) * dx
+ 0.1 * (
theta * u.dx(0)
) * phi.dx(0) * dx
+ ((theta*u)*(theta*u.dx(0))) * phi * dx
)
J = derivative(F, u, du)
solve(
F == 0,
u,
bcs=bc_left,
J=J,
solver_parameters={
"newton_solver": {
"relative_tolerance": 1e-8,
"absolute_tolerance": 1e-8,
"maximum_iterations": 70,
"relaxation_parameter": 0.5,
"linear_solver": "lu"
}
}
)
# Update solution
uold.assign(u)
# -----------------------------------------------
# Store solution after every time step
# -----------------------------------------------
times.append(t)
U.append(u.vector().get_local()[idx])
return x, np.array(times), np.array(U)
# -------------------------------------------------------
if __name__ == "__main__":
Tfinal = 1.0
M = 1500
N = 100
x, times, U = main(
T=Tfinal,
N=N,
M=M,
degree=2
)
print("========================================")
print("Simulation completed successfully.")
print("Number of spatial nodes :", len(x))
print("Number of time levels :", len(times))
print("Time step dt :", Tfinal/M)
print("Solution matrix shape :", U.shape)
print("========================================")
print("Files saved:")
print(" x.txt")
print(" time.txt")
print(" solution.txt")
# ===================================================
# 3D SOLUTION PLOT
# ===================================================
T_mesh, X_mesh = np.meshgrid(
times,
x,
indexing="ij"
)
# ---------------------------------------------------
# Create figure
# ---------------------------------------------------
fig = plt.figure(
figsize=(12, 8)
)
ax = fig.add_subplot(
111,
projection="3d"
)
# ---------------------------------------------------
# Surface plot
# ---------------------------------------------------
surface = ax.plot_surface(
T_mesh,
X_mesh,
U,
cmap="viridis",
edgecolor="none",
antialiased=True
)
# ---------------------------------------------------
# Axis labels
# ---------------------------------------------------
ax.set_xlabel(
"Time $t$",
fontsize=13,
labelpad=10
)
ax.set_ylabel(
"Space $x$",
fontsize=13,
labelpad=10
)
ax.set_zlabel(
"$u(x,t)$",
fontsize=13,
labelpad=10
)
# ---------------------------------------------------
# Title
# ---------------------------------------------------
ax.set_title(
"Finite Element Solution",
fontsize=15,
pad=15
)
# ---------------------------------------------------
# Color bar
# ---------------------------------------------------
fig.colorbar(
surface,
ax=ax,
shrink=0.65,
pad=0.10,
label="$u(x,t)$"
)
# ---------------------------------------------------
# Viewing angle
# ---------------------------------------------------
ax.view_init(
elev=30,
azim=-125
)
plt.tight_layout()
plt.show()
Newton method fails to converge. Any suggestions or examples would be greatly appreciated.