Approximating a maximum value with p-norm function

variational formulation
1

Hello, I’m trying to approximate a maximum value of a field, u_sol, using the p-norm, \int(u_{sol}(x)^pdx)^\frac{1}{p}. I want to be able to take the derivative of the functional so I can eventually use it with optimization (and with stress), which is why I’m using the p-norm as an approximation.

I’ve implemented the p-norm, but I get an unexpected result. I would expect the p-norm of the field to always be larger than the maximum value, but as the p value increases p-norm gets closer to the maximum value. However, I don’t get this behavior from my implementation and I’m not sure what’s going wrong.

So this is a simple example of what I’ve done, I’m using legacy FEniCS 2019.1.0,

from fenics import *
from fenics_adjoint import *
import numpy as np
from ufl import nabla_div

# Define parameters
nelx = 20
nely = 20
nelz = 20
lx = 0.1
ly = 0.1
lz = 0.1
kappa = Constant(20)

# Create a mesh
mesh = BoxMesh.create(
    [Point(0.0, 0.0, 0.0), Point(lx, ly, lz)], # define opposing corners
    [nelx, nely, nelz], # number of elements in each direction
    CellType.Type.hexahedron)
mesh = Mesh(mesh)

# FE function space and functions
V = FunctionSpace(mesh, "CG", 1)
u = TrialFunction(V)
v = TestFunction(V)
u_sol = Function(V)  # function for storing the solution

# Bounday condition
class Left(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0.)
left = Left()
bc = DirichletBC(V, Constant((0.)), left)

# Define forms
f = Constant((1.0781e3))
L = dot(f, v) * dx
a = inner(kappa*grad(u), grad(v)) * dx

# Solve
A, b = assemble_system(a, L, bc)
solver = KrylovSolver("cg", "ilu")
solver.solve(A, u_sol.vector(), b)

# Calculate pnorm
p = 6
J1 = pow(assemble(pow(u_sol, p)*dx(metadata={'quadrature_degree': 5})), 1/p)
p = 20
J2 = pow(assemble(pow(u_sol, p)*dx(metadata={'quadrature_degree': 5})), 1/p)
p = 40
J3 = pow(assemble(pow(u_sol, p)*dx(metadata={'quadrature_degree': 5})), 1/p)

# Compare with pnorms with max value
print('pn-6    =', float(J1))
print('pn-20   =', float(J2))
print('pn-40   =', float(J3))
print('Max val =', u_sol.vector().max())

The output of this is

pn-6    = 0.07120772774142763
pn-20   = 0.17573760154304024
pn-40   = 0.21576352505332388
Max val = 0.2695257644203194

where the first three values are the p-norm evaluated at p=6, p=20, and p=40 and the last value is the maximum value from u_sol. So the p-norm is much smaller than the maximum and it’s not until a very large p value, 40, that it starts to get close.

2

You cannot expect it to converge at a fast rate while limiting the quadrature degree to 5. By increasing the degree, you will approximate the p-th order polynomial better and better (as you are using P-1 as u)

3

Thanks. So I did try increasing the quadrature degree, but it didn’t have that much of an effect on the p-norm when p = 12.

These are the quadrature degrees I tried

# Calculate pnorm
p = 12
J1 = pow(assemble(pow(u_sol, p)*dx(metadata={'quadrature_degree': 5})), 1/p)
J2 = pow(assemble(pow(u_sol, p)*dx(metadata={'quadrature_degree': 10})), 1/p)
J3 = pow(assemble(pow(u_sol, p)*dx(metadata={'quadrature_degree': 12})), 1/p)
J4 = pow(assemble(pow(u_sol, p)*dx(metadata={'quadrature_degree': 15})), 1/p)
J5 = pow(assemble(pow(u_sol, p)*dx(metadata={'quadrature_degree': 20})), 1/p)

# Compare with pnorms with max value
print('\npn-12, q5  =', float(J1))
print('pn-12, q10 =', float(J2))
print('pn-12, q12 =', float(J3))
print('pn-12, q15 =', float(J4))
print('pn-12, q20 =', float(J5))
print('Max val    =', u_sol.vector().max())

And the output

pn-12, q5  = 0.13488788837621457
pn-12, q10 = 0.13488788839202304
pn-12, q15 = 0.134887888392023
pn-12, q20 = 0.134887888392023
pn-12, q40 = 0.13488788839202304
Max val    = 0.2695257644203194

Should I be increasing the quadrature degree even more?

4

Back to your assumption:
Why do you expect p-norm to be larger than the maximum value?
This does not make sense to me, as if you integrate u_max^p, u_max>0 will always be bigger than u^p integrated.