Hello,
I am trying to solve the following fourth-order PDE with Fenics
\nabla^2 \partial_i (z\, \partial_i z) = f
on the unit square. I introduce the following auxiliary fields
\omega_i \equiv \partial_i z
\mu = \partial_i (z \omega_i)
I do this because I cannot differentiate twice fields, whose derivatives are discontinuous across elements. I know that I could have solved this problem with the two fields z and \mu only, but I do need \omega because this is a minimal working example for a more complex problem where \omega exists.
I solve this problem, in a case where I know the exact solution: z = (x^4+y^4)/48, f = (3 x^4 + x^2 y^2 + 3y ^4)/8.
I prepared a minimal working example code for you which reproduces the issue:
from fenics import *
import ufl as ufl
from dolfin import *
import sys
L = 1.0
h = 1.0
function_space_degree = 1
mesh = UnitSquareMesh(32, 32)
i, j, k, l = ufl.indices(4)
xdmffile_err_z = XDMFFile("err_z.xdmf")
xdmffile_err_omega = XDMFFile("err_omega.xdmf")
xdmffile_err_mu = XDMFFile("err_mu.xdmf")
# radius of the smallest cell in the mesh
r_mesh = mesh.hmin()
boundary = 'on_boundary'
n = FacetNormal(mesh)
P_z = FiniteElement('P', triangle, function_space_degree)
P_omega = VectorElement('P', triangle, function_space_degree)
P_mu = FiniteElement('P', triangle, function_space_degree)
element = MixedElement([P_z, P_omega, P_mu])
Q = FunctionSpace(mesh, element)
Q_z = Q.sub(0).collapse()
Q_omega = Q.sub(1).collapse()
Q_mu = Q.sub(2).collapse()
assigner = FunctionAssigner(Q, [Q_z, Q_omega, Q_mu])
class z_exact_expression(UserExpression):
def eval(self, values, x):
# values[0] = np.cos( x[0] + x[1] ) * np.sin( x[0] - x[1] )
values[0] = (x[0] ** 4 + x[1] ** 4) / 48.0
def value_shape(self):
return (1,)
class omega_exact_expression(UserExpression):
def eval(self, values, x):
values[0] = (x[0] ** 3) / 12.0
values[1] = (x[1] ** 3) / 12.0
def value_shape(self):
return (2,)
class mu_exact_expression(UserExpression):
def eval(self, values, x):
values[0] = (7 * x[0] ** 6 + 3 * x[0] ** 4 * x[1] ** 2 + 3 * x[0] ** 2 * x[1] ** 4 + 7 * x[1] ** 6) / 576.0
def value_shape(self):
return (1,)
class f_exact_expression(UserExpression):
def eval(self, values, x):
# values[0] = -16 * (np.cos( 4 * x[0] ) + np.cos( 4 * x[1] ) + np.sin( 2 * x[0] ) * np.sin( 2 * x[1] ))
values[0] = 1 / 8.0 * (3 * x[0] ** 4 + x[0] ** 2 * x[1] ** 2 + 3 * x[1] ** 4)
def value_shape(self):
return (1,)
# Define variational problem
psi = Function(Q)
nu_z, nu_omega, nu_mu = TestFunctions(Q)
z_output = Function(Q_z)
omega_output = Function(Q_omega)
mu_output = Function(Q_mu)
z_exact = Function(Q_z)
omega_exact = Function(Q_omega)
mu_exact = Function(Q_mu)
f = Function(Q_z)
J_Q = TrialFunction(Q)
z, omega, mu = split(psi)
z_exact.interpolate(z_exact_expression(element=Q_z.ufl_element()))
omega_exact.interpolate(omega_exact_expression(element=Q_omega.ufl_element()))
mu_exact.interpolate(mu_exact_expression(element=Q_mu.ufl_element()))
f.interpolate(f_exact_expression(element=Q_z.ufl_element()))
z_profile = Expression('(pow(x[0], 4) + pow(x[1], 4)) / 48.0', element=Q.sub(0).ufl_element())
mu_profile = Expression('(7 * pow(x[0], 6) + 3 * pow(x[0], 4) * pow(x[1], 2) + 3 * pow(x[0], 2) * pow(x[1], 4) + 7 * pow(x[1], 6))/576.0', element=Q.sub(2).ufl_element())
bc_z = DirichletBC(Q.sub(0), z_profile, boundary)
bc_mu = DirichletBC(Q.sub(2), mu_profile, boundary)
assigner.assign(psi, [f, omega_exact, mu_exact])
F_z = ((mu.dx(j)) * (nu_z.dx(j)) + f * nu_z) * dx \
- n[j] * (mu.dx(j)) * nu_z * ds
F_omega = (z * ((nu_omega[i]).dx(i)) + omega[i] * nu_omega[i]) * dx \
- n[i] * z * nu_omega[i] * ds
# F_mu = ((z * omega[i]).dx(i) * nu_mu - mu * nu_mu) * dx
F_mu = (z * omega[i] * (nu_mu.dx(i)) + mu * nu_mu) * dx \
- n[i] * z * omega_exact[i] * nu_mu * ds
F = (F_omega + F_z + F_mu)
bcs = [bc_z]
J = derivative(F, psi, J_Q)
problem = NonlinearVariationalProblem(F, psi, bcs, J)
solver = NonlinearVariationalSolver(problem)
solver.solve()
z_output, omega_output, mu_output = psi.split(deepcopy=True)
xdmffile_err_z.write(project(z_output - z_exact, Q_z), 0)
xdmffile_err_omega.write(project(sqrt((omega_output[i] - omega_exact[i]) * (omega_output[i] - omega_exact[i])), Q_z), 0)
xdmffile_err_mu.write(project(mu_output - mu_exact, Q_mu), 0)
which runs ok
$ python3 solve.py
Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 1.175e+01 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 1: r (abs) = 2.681e-05 (tol = 1.000e-10) r (rel) = 2.282e-06 (tol = 1.000e-09)
Newton iteration 2: r (abs) = 4.876e-07 (tol = 1.000e-10) r (rel) = 4.151e-08 (tol = 1.000e-09)
Newton iteration 3: r (abs) = 1.152e-07 (tol = 1.000e-10) r (rel) = 9.810e-09 (tol = 1.000e-09)
Newton iteration 4: r (abs) = 2.196e-08 (tol = 1.000e-10) r (rel) = 1.869e-09 (tol = 1.000e-09)
Newton iteration 5: r (abs) = 2.963e-09 (tol = 1.000e-10) r (rel) = 2.522e-10 (tol = 1.000e-09)
Newton solver finished in 5 iterations and 5 linear solver iterations.
The Fenics solutions for z and \omega agree with the exact one, but the one for \mu disagrees with the exact one at the corners of the square, see screenshot attached for \mu - \mu_\text{exact}
Why are there such big discrepancies? How may I fix them?
Thank you