Subject: [Ufl] [Question #215432]: True User-defined Operators for UFL

From: "David M. Rogers" <[email-masked]>
Date: Tue, 27 Nov 2012 23:20:59 -0000

New question #215432 on UFL:
https://answers.launchpad.net/ufl/+question/215432

Drs. Alnaes and Logg,

  I would like to create an interface between the Cantera chemical equilibrium code and Fenics in order to solve transport problems in nonideal solutions.  My problem is a nonlinear Poisson equation div ( p(c)  grad c ) = 0.

  The trouble I'm having is in implementing the (horribly complicated, nonlinear) p(c) part, since this involves a call-out to an external code.  The difficulty here is implementing a nonlinear coefficient field, p, that is a function of the object to be solved for rather than simply a function of position.
# nonlin_coef :: (x -> current solution at x) -> (x -> result at x)

 It can be curried, with the type
# nonlin_coef :: (x -> current solution at x) -> x -> result at x
but then will have to be re-compiled for each 'current solution' trial mesh function.

 The correct resolution would be to add a user-defined Expr type that returns a Function object as a function of some other
Function object, but that is currently not possible, since UFL's user guide,  2.12 "User-defined operators" can only be defined
with respect to already existing UFL operators.

  This leads to the advice in UFL's user guide sec. 4.6 "Extending UFL", which appears less hopeful.  In theory, it should be possible to implement this sort of functionality by requiring a series of callback functions from the user to represent successive derivatives.  How difficult would this be to implement?

  Also, the networkx toolkit (http://networkx.lanl.gov/examples/drawing/unix_email.html) has some really nice options for displaying graphs.

Sincerely,
~ David M. Rogers
  Postdoc, Nanobiology Department
  Sandia National Labs

-- 
You received this question notification because you are a member of UFL
Team, which is an answer contact for UFL.

---

From: Martin Sandve Alnæs<[email-masked]>
Date: Wed, 28 Nov 2012 09:06:07 -0000

Question #215432 on UFL changed:
https://answers.launchpad.net/ufl/+question/215432

    Status: Open => Answered

Martin Sandve Alnæs proposed the following answer:
I'm not quite sure what you're asking for here, but I'll try to point
you in the right direction. Are you aware of the distinction between UFL
and DOLFIN? You will have to evaluate your external code and its
derivatives behind the appropriate DOLFIN interfaces, and use
placeholder functions for the results in UFL.

The tricky part is not evaluate the nonlinear coefficient field, that
you can just implement in C++ behind dolfin interfaces, but the
linearization. What is your linearization strategy? You don't specify
what kind of derivatives you can get for p(c). If you can compute

s: Omega -> R
s = \frac{\partial p(c)}{\partial c}

then you can implement the derivative of your weak form as

c = Coefficient(V)
s = Coefficient(V) # represents s = dp/dc
q = Coefficient(V) # represents q = p(c)
ct = TrialFunction(V)
v = TestFunction(V)
a = s*ct*dot(grad(c), grad(v))*dx + q*dot(grad(ct), grad(v))*dx

(not verified, so don't trust this!)

To use the Newton solver in dolfin, you need to implement the NonlinearProblem as in this example:
http://fenicsproject.org/documentation/dolfin/1.0.0/cpp/demo/pde/cahn-hilliard/cpp/documentation.html
and remember to update q = p(c) and s = dp/dc before each assemble.

---

From: "David M. Rogers" <[email-masked]>
Date: Wed, 28 Nov 2012 19:50:56 -0000

Status: Answered => Open

David M. Rogers is still having a problem:
  Thanks for the quick reply.  I was thinking along those lines (and
would of course need dp/dc), but was unsure where to start.  Are you
saying that I can do the following to the CahnHilliard example?

c = Function(V) # to be solved for, index J
p = Coefficient(V) # represents q = p(c), indices I,J
v = TestFunction(V), index I

L = dot(p*nabla_grad(c), nabla_grad(v))*dx

ct = TrialFunction(V) # index K
dpdc = Coefficient(V) # represents s = dp/dc, indices I,J,K

# a = derivative(L, c, ct) # Gateaux derivative in dir. of ct
# From Sec. 1.2.4 of the Fenics Manual:
a = dot(p*nabla_grad(ct), nabla_grad(v))*dx \
  + dot(dpdc*ct*nabla_grad(x), nabla_grad(v))*dx

class UserFunction(Expression): ...
class UserDeriv(Expression): ...

fp = UserFunction(c) # stores location of c's data
fdpdc = UserDeriv(c) # stores location of c's data

# Class for interfacing with the Newton solver
class CahnHilliardEquation(NonlinearProblem):
    def __init__(self, a, L, p, dpdc):
        NonlinearProblem.__init__(self)
        self.L = L
        self.a = a
        self.p = p
        self.dpdc = dpdc
        self.reset_sparsity = True
    def F(self, b, x):
        self.p.interpolate(f) # sees new 'c' when called?
        assemble(self.L, tensor=b)
    def J(self, A, x):
        self.p.interpolate(fp)
        self.dpdc.interpolate(fdpdc)
        assemble(self.a, tensor=A, reset_sparsity=self.reset_sparsity)
        self.reset_sparsity = False

  That would probably make the best solution as long as JIT compiling
will not be re-done.  Since  I'm defining p(c) and dpdc as 'Coefficient'
fields, this is a Dolfin issue, so I'll send further questions their
way.  If it works I'll try and send the example back to the Fenics
project.

Thanks!
~ David

---

From: Launchpad Janitor <[email-masked]>
Date: Fri, 14 Dec 2012 19:51:43 -0000

Status: Open => Expired

Launchpad Janitor expired the question:
This question was expired because it remained in the 'Open' state
without activity for the last 15 days.