How to write nonlinear function from data in terms of solving variable in nonlinear variational problem

Question

I'm trying to solve a nonlinear problem, say
nabla.(k(u).nabla(u)) = C.

I have a problem in that i don't know k(u) in terms of a function of u, or a function of spatial coordinates, but only as a list of u and its corresponding k. I currently have a function which interpolates the data, and if a solution u is known, I can find out the corresponding k's from this function.

However, I don't think this is working once it's incorporated into the problem as it solves. I've been told making a subclass in Expression might be worth looking into, but i have no idea how i'm supposed to do this since i don't know things in terms of x[0] etc.

Any help appreciated - below is a snippet of code which explains where my k is.
interpolate_table is my interpolating code that takes points on u and uses k_data which is teh list of data and outputs corresponding k (if i do plot(project(k,V)) it works as expected). I had used this code with k_data being points with k=1+u^2 and C was the corresponding source term with u=x^2 as the solution.

u = Function(V)
v = TestFunction(V)

k = interpolate_table(u,k_data)


C = Expression("-2.0-10.0*x[0]*x[0]*x[0]*x[0]")
diff_part = inner(k*grad(u),grad(v))*dx
source_part = S*v*dx

Res = diff_part - source_part
J = derivative(Res,u)
problem = NonlinearVariationalProblem(Res, u, bcs, J)
solver  = NonlinearVariationalSolver(problem)
solver.solve()

prm = solver.parameters['newton_solver']
prm['relative_tolerance'] = 1E-6
prm['maximum_iterations'] = 100

Answer 1

How about replacing k = interpolate_table(u,k_data) by the following?

class K(Expression):
    def eval(self, value, x):
        value[0] = interpolate_table(u(x), k_data)
k = K()

I don't think FEniCS knows about the dependence of k on u, so instead of using J = derivative(Res, u) you might have to derive and code the expression for the Jacobian yourself.

Comments

Hi, thanks for this.

I tried this and whilst it works for 1D, it gets stuck at the second timestep in 2D. it is also very slow...

the 2d solver gets stuck at the value[0] = interpolate_table(u(x),k_data) bit. Do you have any ideas why this would be? Since i don't know k in terms of spatial dimensions, i dont' know how else to write the expression...

Many thanks for your help

---

It is slow since K() is a Python expression that evaluates slower than its C counterpart. The following question might be of help, although it is only partly answered:
http://fenicsproject.org/qa/8227/jit-expression?show=8227#q8227

I don't see what is wrong with the 2D case. You will more likely get a reply if you post a minimal code example that produces the error.