Conditioning Finite-element solution to Non-negative values

1

Hello, FEniCS community,

I am stuck with what seems to be a very simple problem. I have a dynamic simulation (with 6 PDEs), and for some numerical reasons, I am obtaining negative values within my solution (u). I need to force the solution to be non-negative. That is, if there is any negative value within the spatial function, it must be transformed into zero. For this purpose, I am using UFL conditionals as shown below. However, I get the following error:

Expecting scalar arguments.
Traceback (most recent call last):
File “Variational_Problem_EF.py”, line 284, in SafeU
return conditional(gt(x, DOLFIN_EPS), x, DOLFIN_EPS)
ufl.log.UFLException: Expecting scalar arguments.

So I also tried to split u into each of the six Finite element solutions suchlike:

u_1, u_2, u_3, u_4, u_5, u_6 = split(u)
conditional(gt(u_1, DOLFIN_EPS), u_1, DOLFIN_EPS)
conditional(gt(u_2, DOLFIN_EPS), u_2, DOLFIN_EPS)
conditional(gt(u_3, DOLFIN_EPS), u_3, DOLFIN_EPS)
...

However, it seems that each of these lines of code create a C++ expression. Further, these expressions must be compiled to actually update the solution by eliminating all negative values to zeros.

I post an exemplification of the code, notwithstanding, it cannot be executed, it is only aimed to facilitate the explanation above.

P1 = FiniteElement('P', triangle, 1)
element = MixedElement([P1, P1, P1, P1, P1, P1])
V = FunctionSpace(mesh, element)

# Define test functions
v_1, v_2, v_3, v_4, v_5, v_6 = TestFunctions(V)
# Define functions for concentrations and temperature
u = Function(V)
u_n = Function(V)

# Split system functions to access components
u_n1, u_n2, u_n3, u_n4, u_n5, u_n6 = split(u_n)
u_1, u_2, u_3, u_4, u_5, u_6 = split(u)

# Variational formulation
F = ...

for n in range(num_steps):
    t += deltaT  # Update current time

    solve(F == 0, u, bcs, solver_parameters={"newton_solver": {"relative_tolerance": 1e-5}})

    def Non_zeroU(x):
        return conditional(gt(x, DOLFIN_EPS), x, DOLFIN_EPS)

    u_n.assign(Non_zeroU(u))

Thank you so much for any suggestions or comments on this topic.

All the best,
Santiago

2

ufl.conditional is a computational symbolic expression of the condition provided. It will not apply a numerical operator to a vector of degrees of freedom as your call to u_n.assign is attempting.

There are two approaches I would investigate:

  1. Your problem is already set up as if it were nonlinear: try solving it with u_1 through u_6 replaced by your conditional symbolic representations. This will probably be difficult (maybe impossible) for the Newton solver to resolve without an appropriate initial guess though.
  2. Get the local values of the underlying vector of u, and manipulate them as you wish between time steps. E.g. do_something_with(as_backend_type(u.vector()).get_local()).
1 Like
3

Hello nate

Thank you for the fast reply.

Until now I have tried the second approach:

for n in range(num_steps):
    t += deltaT  # Update current time
    solve(F == 0, u, bcs, solver_parameters={"newton_solver": {"relative_tolerance": 1e-5}})
    # Manipulating data in solution
    Uvector = as_backend_type(u_n.vector()).get_local()
    Uvector[Uvector <= DOLFIN_EPS] = DOLFIN_EPS
    # Reassigning fixed solution
    u_n.assign(...)

It is possible to manipulate data as wish, but now I need to reassign Uvector to u_n.
Any idea about how this could be done?

Best regards,
Santiago

4

Hello, you should do something like

u_n.vector().set_local(Uvector)
5

@nate what is the syntax for achieving non-negativity in dolfinx? I have a solution from ("CG", 1) that I want to ensure is positive semi-definite.

6

u.x.array[u.x.array<0] = 0 for instance

1 Like