Non-Linear Coeffecients: dependent on trial functions;

variational formulation
1

Hi All,

I am working on PDE of the below weak form: (not the exact problem I am dealing with, but I have simplified it without worrying about its physical sense)

a = (p*S(p)-U(p))*inner(epsilon(v),Identity(ndim))*dx + q*tr(epsilon(u))*dx

u,p are trial functions
v,q are test functions

I am handling this problem using a mixed finite element. (u is vector and p is scalar solutions) Correspondingly the mixed element is P2-P1.

Now in the first term of the PDE the coefficients S(p) and U(p) are making the problem nonlinear.

  1. Nevertheless, I have defined the function S(pf) as below:

    def S(pf):
    return ((-pf)^(1/2)+1)^-1

However, I am unable to do this and the error is as follows:

UFLException: Cannot take the power of a non-scalar expression <ComponentTensor id=140377318371880>

  1. U(pf) is not so straight forward. It is a definite integral of the function S(pf) in the limits p and 0. I am not sure how to handle this. I have tried using sympy as below:

    def U(pf):
    z = sympy.Symbol("z")
    return sympy.integrate(S(z),(z,0,pf)) - pf*S(pf)

With this I get an error as below when evaluating the return of the U(pf) function.
NotImplementedError: Cannot take length of non-vector expression.

I have tried another way using the assemble function in FEniCS as below:

    def U(pf):
        mesh = IntervalMesh(10,pf,0.0)
        z = SpatialCoordinate(mesh)
        return assemble(S(z)*dz(degree=2)) - pf*S(pf)

With the second method I get an error on the evaluation of the IntervalMesh definition as below:

TypeError: __init__(): incompatible constructor arguments. The following argument types are supported:
    1. dolfin.cpp.generation.IntervalMesh(arg0: int, arg1: float, arg2: float)
    2. dolfin.cpp.generation.IntervalMesh(arg0: MPICommWrapper, arg1: int, arg2: float, arg3: float)

Invoked with: 10, Indexed(Argument(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 3), MixedElement(VectorElement(FiniteElement('Lagrange', triangle, 2), dim=2), FiniteElement('Lagrange', triangle, 1), FiniteElement('Lagrange', triangle, 1), FiniteElement('Lagrange', triangle, 1))), 1, None), MultiIndex((FixedIndex(2),))), 0.0

I believe these errors are related to the solution p being passed directly for function evaluation.
The functions S(pf) and U(pf) need to be evaluated at every node i.e. for every value of p in the domain and should be returned to the weak form.

I believe that I am making a basic mistake and cannot get my head around it. Can anyone of you help with this?

Many Thanks,
Sid

2

I have tried today another way to evaluate the coefficient function U(pf) that is dependent on the trial function p.

I tried to ensure that the return value into the weak form is of the form UFL. In the coefficient function definition I looped on the vertex values of the trial function p (using p.compute_vertex_values()) and within each loop I evaluated the expression that I need (integration of a function in my case) and then assigned the computed value to a corresponding vertex value of the coefficient function which was apriori defined as a function(maybe I should try as a trial function) of the mesh domain similar to the actual trial function p. This seems to be logical and intuitive but doesn’t seem to work as such for a non-linear problem. However, this kind of function evaluation by looping seems to work if loop over the vertex values of a previous solution (i.e. a function of instead of a trial function like p). So end of the day I probably have to resort to some kind of iterative approach (Picard or Newton) or approximate that my coefficients need to be evaluated from a previous know solution. Nevertheless, please see what I have done and the corresponding error in the case of direct non-linear approach.

Q       = FunctionSpace(mesh,'P',1)
Uf  = Function(Q)
def U(pf):
    coordinates = mesh.coordinates()
    vertex_values_p = pf.compute_vertex_values()
    vertex_values_U = Uf.compute_vertex_values()
    for i,x in enumerate(coordinates):
        z = sympy.Symbol("z")
        vertex_values_U[i] = sympy.integrate((((-z)^(1/2)+1)^-1),(z,0,vertex_values_p[i])) - vertex_values_p[i]*(((-vertex_values_p[i]/M)**(1/pm)+1)**-1)
    return Uf

The error is as below when evaluating the function call:

215     coordinates = mesh.coordinates()
--> 216     vertex_values_p = pf.compute_vertex_values()
217     vertex_values_U = Uf.compute_vertex_values()

AttributeError: 'Indexed' object has no attribute 'compute_vertex_values'
fenics@2bc3e4917bf0:~/shared$

Basically it says that the when I pass the trial function as an argument (nonlinear case) it says that there exist no attribute called vertex values but when I do the same by passing a function (previous solution in linear case) it extracts nicely the vertex values.

The weird part is in my original post the error due to Coefficient function S(pf) was because I recalled it in the U(pf). As such S(pf) does not seem to have any problem even if I try to evaluate it in a non linear case using trial function p. Only U(pf) seems to have this issue. I understood this when I stopped calling S(pf) in U(pf) and instead directly used the function expression of S(pf) and additional thing I did was a definite integration. I believe the problem lies in this additional integration.

Any ideas would be appreciated.

many Thanks,
Sid

3

Hey Siddhartha, I found your query from your post here.

So p is a trial function and is a UFL object. Your definition of S(p) would need UFL operators to be valid.

def S(pf):
   return pow((sqrt(-pf)+1),-1)

(FYI : pow and sqrt are from dolfin or ufl namespaces and not python math ones)

As for U(p), since it is a definite integral, I would consider it a UFL Constant which is a function of a trial function p. I am not too sure what you are trying to implement here - nor the skills to solve what you need. Just useful info I can give you : It should be possible to use assemble() over a form to obtain the definite integral.

Some Ref for valid UFL forms and Operators can be found in documentation here : https://fenics.readthedocs.io/projects/ufl/en/latest/

You can maybe post this at UFL mailing list for better ideas from experts

Cheers.