improving precision of conditional Max function (truncation)

Question

I'd like to flag when a function has exceeded a threshold, equivalent to truncating a function.
For instance I can apply Max(f,Constant(0)), which does the job of setting all negative values to 0.

The problem is that the precision of that operation is not necessarily high (unless the mesh is particularly fine). I illustrate the problem in the script below, truncating the negative values of cos(x). The result of the Max operation is quite noisy, and even returns some negative values. Ideally I want it to simply pass the value of the function of the original function where it is positive, and strictly zero otherwise.

Increasing the number of points in the mesh does help, but I guess it's not always convenient to do that in production runs. Is there a better, cleaner method of truncating function values?

from dolfin import *

mesh = IntervalMesh(20,0,10)
V = FunctionSpace(mesh, "CG", 1)

fe = Expression("cos(x[0])")
ff = interpolate(fe,V)

truncated_f = project(Max(ff,Constant(0)),V)

print("original: ", ff.vector().array())
print("truncated: ", truncated_f.vector().array())

plot(ff,title="Original")
plot(truncated_f,interactive=True,title="Truncated")

Comments

What you perceive as an inaccuracy in the Max operator is actually caused by the projection.

---

Interesting point. That means I don't need to worry about the "noise" during the actual solving of the PDE, doesn't it? The noise will only show up at the end when I try to access the values of the function coming from the final solution.

---

Projection is an operation that in general does not preserve pointwise properties.

Think of max(f, 0) as a function taking pointwise non-negative values. In general, max(f,0) will not be in V, even if f is. With the projection, you find its best approximation (in L^2 metric) in V. This approximation is not in general non-negative.

---

I think that makes sense. Essentially the projection is attempting to approximate a step function, with resultant "ringing" at the edge of the step.

Answer 1

You can do the following:

truncated_f = Function(ff)
truncated_f.vector()[ff.vector()<0] = 0.0

Comments

Thanks Øyvind.

That solution works when ff is an ordinary Function.

But it doesn't work when setting up a weak formulation that contains a truncation of the solution, i.e. when ff is a TrialFunction (or a form containing the trial function). In this case Function(ff) gives an error like

  File "/usr/lib/python2.7/dist-packages/dolfin/functions/function.py", line 287, in __init__
    raise TypeError("expected a FunctionSpace or a Function as argument 1")
TypeError: expected a FunctionSpace or a Function as argument 1

The same error happens if ff is a form, e.g. injecting ff = sin(ff) into the test script. The example script will work if I then also project that form (ff=project(ff,V)) before taking your truncated_f = Function(ff). But again, the projection operator doesn't work if ff is a TrialFunction.