Maximum field values: the "right" approach

Question

Dear all,

I'd like to retrieve the maximum field value (or a function of it) over my domain. Right now I'm using the suggested way of projecting and using max:

solve(A, u.vector(), b, 'lu')

def u_norm(v):
    return sqrt(v[0] ** 2 + v[1] ** 2)
print norm_u.vector().max()

I don't know if this is guaranteed to result in the "real" maximum value of u. Can I assume safely that this is the correct approach?

Cheers!

Answer 1

Why not just use

norm(u, 'linf') 

? Isn't that the same thing (but designed for MPI considerations too).
I think you can also use

u.vector().norm('linf') 

(untested)
but isn't MPI'd (you get the maximum value on each node, not the maximum value overall).

Comments

Good point, for a global max you need a call MPI::max to collect the max value.

---

Awesome! I was looking also for MPI solutions, thanks! :)

So how can I retrieve the maximum norm of a field? The norm gives me the norm of u, not vector of norms of u, as I understand.

---

To get a global max/min, you can use the dolfin helper functions

double max_u = u->vector()->max();   // get the local max
max_u = dolfin::MPI::max(MPI_COMM_WORLD, max_u);  // get the global max

---

Ah, okay, now I think I understand what you want. You have a vector field and you want to know the maximum magnitude right?

Firstly, norm(u,'linf') is incorrect. Actually it should be norm(u.vector(),'linf') (the former didn't work on my computer). What it is doing in that case is looking at the solution vector (without any knowledge that that solution vector in fact represents a vector field, not 2 scalar fields) and returning the linf norm (the maximum value). In effect, you are finding the maximum component, not the maximum magnitude.

So then, from what I can think of, your original method is pretty much as good as it's going to get. viz: you have to get a vector of the magnitudes and then take the linf norm of that.

BTW: Instead of your u_norm function , you can also just do

inner(u,u)

although I don't know if it would be any faster or anything.

Answer 2

It will be the maximum value at a DOF. The vector, norm_u.vector(), is a vector of values for each DOF, from that vector the maximum value is taken. This isn't the same as the maximum value of the Function itself unless you have linear elements.

If you want to measure "accurately" the maximum value of the vector u, you might want to interpolate onto a finer mesh with linear elements, and take the maximum value from there. This might be a convenient way to improve that resolution since the interpolation is fairly fast, and the search for a max value is also fast.

If you really need it to be accurate, you might want to use a root finder of some type to find the location of the maximum value and then evaluate that point, which will consider the basis functions. This approach requires a lot of point evaluations and of course a root solver.

Comments

Thanks for the pointers, Charles! So with an interpolation, in Fenics it would be like

mesh = RectangleMesh(0, 0, Width, Height, 10, 20)

# Solve u on mesh...

finerMesh = RectangleMesh(0, 0, Width, Height, 100, 200)

V_2D = TensorFunctionSpace(finerMesh, "CG", 1)
new_u = project(u_norm(u), V = V_2D)

print '> Max norm of displacement: ', new_u.vector().max()

What about the root finder? Is some implementation available in Fenics?

---

With me not generally working in Python, that looks right.

For the root finder, I don't think there is anything built in (thought, I haven't looked for it specifically). There are a lot of approaches you could use to do it... probably no best answer (global/local searches for extrema). I like the idea of the root finder as it abstracts away details of the function space. This might be a tangent not worth pursuing.

In general, if you know something of the basis functions themselves, you should be able to calculate a max value in the element knowing the DOFs. Taking a quick look there are methods to calculate the basis function derivatives - which perhaps you could use to find where in each cell to check for a maximum?

What you have above is probably good enough, considering the function space is not exact anyhow?