Different results when using interpolate, project and expression

Hi, I have met a strange problem when solving time dependent non-linear Poisson form equation with the Picard iteration. The equation includes several highly non-linear parameters. Here I use k[f(u)] to represent the parameter, where u is the unknown variable, f(u) is a function of u, and k is a function of f(u). As the Picard iteration refers to assigning a already known uk into the parameter k and get k[f(uk)] in each iteration, I’ve tried three ways to define and update the k[f(uk)]:

  1. define and update f and k as Expression:
    f_ex=Expression('...u.....',degree=2,u=uk)
    k_ex=Expression('...f.....',degree=2,f=f_ex)
    k=interpolate(k_ex,V)
    while ....:
        uk.assign(ux) #ux refers the solution
        f_ex.u=uk
        k_ex.f=f_ex
        k.assign(interpolate(k_ex,V))

This method leads to satisfying stable iteration.

  1. define and update f and k as:
    f_ex=Expression('...u.....',degree=2,u=uk)
    k_ex=Expression('...f.....',degree=2,f=f_ex)
    k=interpolate(k_ex,V)
    while ....:
        uk.assign(ux)
        f_ex.u=uk
        k_ex.f=interpolate(f_ex,V)
        k.assign(interpolate(k_ex,V))

The only difference between this method and the method 1 is the interpolated f is used to get k_ex but not the f_ex. Howver, this leads to finally unstable iterations (nan results).

  1. define and update f and k by project:
    k=project(f(uk),V)
    while ....:
        uk.assign(ux)
        k.assign(project(f(uk),V))

However, this also leads to finally unstable iterations (nan results).
Why these methods which shoud have same k but lead to different iteration results even with nan results while method 1 iterated successfully?

1 Like

This is highly dependent on your choice of function space. As you are using the degree=2 argument in the expression, this will be interpolated into a second order function space. However, if your space V is a first order space, your second interpolation will reduce the accuracy of f.
Similarly, projection and interpolation does to different kinds of approximation: interpolation sets values exactly at degrees of freedom, minimizes the integrated difference over the mesh.

For me to be of any further help, you would have create a minimal code example, as described in: Read before posting: How do I get my question answered?

2 Likes