Ah, I see! (you need a + not a - by the way).

That would come from a DG discretisation, or finite volume at lowest order. The only difference between DG and finite volume is that in FV the volume integral vanishes, whilst the surface flux integrals remain (which would be constructed in the way that David described in Firedrake).

all the best

--cjc

On 3 October 2015 at 09:23, Eike Mueller <E.Mueller@bath.ac.uk> wrote:

Hi Colin,

it's coming from the implicit discretisation of the mass conservation equation:

d\rho/dt - \div(\rho*u) = 0.

Thanks,

Eike

Sent from my iPad

On 3 Oct 2015, at 09:07, Colin Cotter <colin.cotter@imperial.ac.uk> wrote:

Hi Eike,
  It still very much depends on what equation is being solved in what space, and how it has been linearised. If you tell us which part of the model it is from then I can tell you how we do it in our codes.
All the best
Cjc

On 3 Oct 2015 08:09, "Eike Mueller" <E.Mueller@bath.ac.uk> wrote:

Hi David,

I'm trying to understand some code in DYNAMO, where it is done as I described (at least that's what I think - I'm reverse engineering the existing code, and waiting for clarification from the Met Office). I was hoping to shed some light on this by finding out how fluxes are done in firedrake.

Ultimately, I want to assemble the weak form of the linear operator 

u -> A(u) = \div F with F=\phi*u

i.e. <\psi, \div(\phi*u)> with u \in Hdiv, \psi, \phi \in L2. \psi \in L2 is a test function. \phi is a background reference field, and in a sense it doesn't really matter how the flux is constructed, as long as \phi*u is a good approximation to the continuum expression (this is because the linear operator shows up in the linearisation of the non-linear problem). The action of A is implemented by constructing F as described and then looping over all cells and calculating the divergence as <\theta,\div F>. Would a more natural definition of the flux be to find F such that

phi*dot(u,v)*dx = dot(F,v)*dx

for all test for functions v? The thing is, I want to really assemble a linear operator A, so having a solve in there wouldn't work.

Thanks,

Eike

Sent from my iPad

On 3 Oct 2015, at 04:36, David Ham <David.Ham@imperial.ac.uk> wrote:

Hi Eike,

Can you give some context. What you have described is a very unfiniteelement concept. The usual FE concept of the flux is the integral over a facet:

phi*dot(u, n) * ds

Where n is the FacetNormal. In that case everything gets evaluated at quadratic equation points so the issue you raised does not occur.

Cheers,

David

On Fri, 2 Oct 2015 at 20:21 Eike Mueller <E.Mueller@bath.ac.uk> wrote:

Dear firedrakers

if I have a scalar field phi and a vector field u, how can I calculate the flux f=phi*u in firedrake? I need to be able to loop over the dofs i of u, and for each one calculate f_i = sum_k phi_k*e_k(x_i)*u_i where e_k(x_i) is the k-th basis function of phi evaluated at the node x_i of u. But I somehow have to tell it that I want phi evaluated at u's nodes, not the other way round.

Thanks,

Eike

Sent from my iPad
_______________________________________________
firedrake mailing list
firedrake@imperial.ac.uk
https://mailman.ic.ac.uk/mailman/listinfo/firedrake

_______________________________________________
firedrake mailing list
firedrake@imperial.ac.uk
https://mailman.ic.ac.uk/mailman/listinfo/firedrake

_______________________________________________
firedrake mailing list
firedrake@imperial.ac.uk
https://mailman.ic.ac.uk/mailman/listinfo/firedrake

--

http://www.imperial.ac.uk/people/colin.cotter

www.cambridge.org/9781107663916