First order derivatives in non-linear problem

dolfinx
, ,
1

Dear FenicsX users,

I am trying to solve the following non-linear problem:

image

On the [0,1] one-dimensional interval with homogeneous Neumann boundary conditions for u.

I have been able to successfully solve other non-linear problems, but this one, despite being quite simple, yields the following error:

image
image1459×180 9.25 KB

And after playing around a bit, I found out the problem is the first order derivative term in the variational formulation. With “u” being the name of my solution and “v” that of the test function, I am using “ufl.inner(u.dx(0),v)” (full code below). How may I get around the error? Thanks in advance!

import os
import numpy as np
import ufl
from dolfinx import log, plot
from dolfinx.fem import Function, FunctionSpace
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.io import XDMFFile
from dolfinx.mesh import CellType, create_unit_square, create_interval
from dolfinx.nls.petsc import NewtonSolver
from ufl import dx, grad, inner
from dolfinx import fem, mesh
from mpi4py import MPI
from petsc4py import PETSc

import matplotlib.pyplot as plt
import scipy

comm = MPI.COMM_WORLD
N = 1000
r_min = 0
r_max = 1
msh = create_interval(comm, N, [r_min, r_max])
deg = 1
V = FunctionSpace(msh, ("Lagrange", deg))

# Trial and test functions
v = ufl.TestFunction(V)

u = Function(V)  # solution

# Non linearity
fun = u**2

SRC = Function(V)
SRC.interpolate(lambda x: 1 + x[0] - x[0])

# Weak statement of the equation

F = - inner(u.dx(0),v) * dx - inner(SRC, v) * dx + inner(fun,v)*dx

problem = NonlinearProblem(F, u)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = 1e-6

# We can customize the linear solver used inside the NewtonSolver by
# modifying the PETSc options
ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "preonly"
opts[f"{option_prefix}pc_type"] = "lu"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()

r_plot = np.linspace(r_min,r_max,deg*(N +1) )

r = solver.solve(u)
print(f"num iterations: {r[0]}")

#Plot
uvalues = u.x.array
plt.plot(r_plot, uvalues, label='Solution')
plt.grid(True)
plt.tight_layout()
plt.show()```