Re: [firedrake] Flux calculation
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
Is Dynamo currently using a fixed reference state or is it using the state from the previous tilmestep for linearisation? We are currently ducking some of these issues by using a fixed reference state. cheers --cjc On 3 October 2015 at 09:44, Colin Cotter <colin.cotter@imperial.ac.uk> wrote:
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
-- http://www.imperial.ac.uk/people/colin.cotter www.cambridge.org/9781107663916
Hi Colin, thanks. Yes, I can see how one would want to use either DG or FV for this, since the field \rho is discontinuous, and hence \rho*u is also discontinuous. If I assume for a moment that \rho = 1, then I should just be able to just construct the weak form of the operator, which is <\psi,\div(u)>, and this becomes psi*div(u)*dx (with \psi being a test function). I think what is happening in the code is that for \rho*u an approximation F \in Hdiv is constructed (this requires the evaluation of \rho at the nodes of u), and then this is used to construct psi*div(F)*dx via an integral over the cell. Does that make sense? The Met Office uses the previous state, as far as I can tell. How would a fixed reference state help? Is it because you can then express \rho as a function in H1? Eike Sent from my iPad On 3 Oct 2015, at 09:44, Colin Cotter <colin.cotter@imperial.ac.uk<mailto:colin.cotter@imperial.ac.uk>> wrote: 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<mailto: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<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 -- http://www.imperial.ac.uk/people/colin.cotter www.cambridge.org/9781107663916<http://www.cambridge.org/9781107663916> [http://assets.cambridge.org/97811076/63916/cover/9781107663916.jpg] _______________________________________________ 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