Re: [firedrake] Flux calculation
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
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<mailto: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<mailto: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<mailto: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<mailto: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<mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake _______________________________________________ firedrake mailing list firedrake@imperial.ac.uk<mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake _______________________________________________ firedrake mailing list firedrake@imperial.ac.uk<mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake
participants (2)
-
Colin Cotter
-
Eike Mueller