Subject: [Ufl] [Question #219998]: unconventional linear and bilinear forms

From: Paul Constantine <[email-masked]>
Date: Thu, 24 Jan 2013 04:21:00 -0000

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

I'm wondering if the following is possible to express easily in UFL. I'll write it in math/latex style. 

Let's say my domain is an interval x \in [a,b]. Now let u = u(x) \in R^n and v = v(x) \in R^n. These are n-vectors whose elements depend on x. And let A = A(x) \in R^n x R^n be a matrix whose elements depend on x. 

I would like to work with linear forms such as

\int_a^b u^T v dx

and bilinear forms such as

\int_a^b u^T A v dx.

The matrix A might be a "mesh function" that returns an nxn matrix for each mesh point.

Reading through the FEniCS book -- particularly Chapter 17 -- I can't tell if such a thing is possible.

-- 
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: Thu, 24 Jan 2013 06:11:23 -0000

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

    Status: Open => Answered

Martin Sandve Alnæs proposed the following answer:
What is n?

Martin
Den 24. jan. 2013 05:26 skrev "Paul Constantine" <
[email-masked]> følgende:

> Question #219998 on UFL changed:
> Paul Constantine gave more information on the question:
> Oh. I also want to be able to take du/dx \in R^n.
>
> --
> _______________________________________________
> Mailing list: https://launchpad.net/~ufl
> Post to     : [email-masked]
> Unsubscribe : https://launchpad.net/~ufl
> More help   : https://help.launchpad.net/ListHelp
>

---

From: Marie Rognes <[email-masked]>
Date: Thu, 24 Jan 2013 09:11:14 -0000

Status: Open => Answered

Marie Rognes proposed the following answer:

Yes, this is perfectly feasible and  not so unconventional (at least for moderate n): Here's some suggested DOLFIN code (replace VectorFunctionSpace with VectorElement etc to get UFL) 

from dolfin import *
mesh = UnitIntervalMesh(20)

n = 30
V = VectorFunctionSpace(mesh, "CG", 1, dim=n)
u = TrialFunction(V)
v = TestFunction(V)

a = inner(u, v)*dx

Q = FunctionSpace(mesh, "DG", 0)
f = Function(Q)
A = as_matrix([[f for i in range(n)] for i in range(n)])

a = inner(A*u, v)*dx

---

From: Martin Sandve Alnæs<[email-masked]>
Date: Thu, 24 Jan 2013 09:21:03 -0000

Martin Sandve Alnæs proposed the following answer:
As Marie says, but for 3 billion forget it and find another approach. Is
that even possible? You have a O(n^2) computation there...

Martin
Den 24. jan. 2013 07:51 skrev "Paul Constantine" <
[email-masked]> følgende:

>     Status: Answered => Open
> Paul Constantine is still having a problem:
> I have systems with n=3 and n=3 billion. I'll start wherever.
>
> --

---

From: Paul Constantine <[email-masked]>
Date: Thu, 24 Jan 2013 16:15:51 -0000

Paul Constantine posted a new comment:
Thanks! That's tremendously helpful. This comes from a recent idea for
parallel-in-time simulations of chaotic dynamical systems. See, e.g.,
http://arxiv.org/abs/1211.2437

In this case n=3 would be something like the Lorenz system, while n=3e9
would be the spatially discretized, linearized forward and adjoint
operators from a Navier-Stokes model at high Reynolds number (turbulent
flow). Obviously, the latter is out of reach in practice. But that's
what we put on proposals. :) However, something like the Kuramoto-
Sivashinsky equation would be the next step after Lorenz. In this case,
the matrix A would still be a spatially discretized operator, but its
physical dimension would be 1 and n=100 would be a reasonable
discretization.

Maybe there's a way to write it in UFL that would let DOLFIN do the
spatial discretization automatically. Can you express one partial
derivative for a v=Function(V) where V=FunctionSpace(<2d mesh>)?

---

From: Jan Blechta <[email-masked]>
Date: Sun, 24 Feb 2013 16:01:10 -0000

Jan Blechta proposed the following answer:
> Can you express one partial derivative for a v=Function(V)
> where V=FunctionSpace(<2d mesh>)?

What about v.dx(i) ?

---

From: Paul Constantine <[email-masked]>
Date: Sun, 24 Feb 2013 16:06:00 -0000

Status: Answered => Solved

Paul Constantine confirmed that the question is solved:
Thanks Jan Blechta, that solved my question.