- firedrake - imperial.ac.uk

Re: [firedrake] Stochastic differential equations
by Colin Cotter 28 Jun '16

28 Jun '16

Hi Francis, The model is described in this paper: https://arxiv.org/abs/1410.8311 all the best --cjc On 27 June 2016 at 14:45, Francis Poulin <fpoulin(a)uwaterloo.ca> wrote: > Hello Colin, > > I would be curious to hear more about the stochastic QG model. Should I > ask Ali about this? > > Cheers, Francis > > > > ------------------ > Francis Poulin > Associate Professor > Department of Applied Mathematics > University of Waterloo > > email: fpoulin(a)uwaterloo.ca > Web: https://uwaterloo.ca/poulin-research-group/ > Telephone: +1 519 888 4567 x32637 > > ------------------------------ > *From:* firedrake-bounces(a)imperial.ac.uk [firedrake-bounces(a)imperial.ac.uk] > on behalf of Colin Cotter [colin.cotter(a)imperial.ac.uk] > *Sent:* Monday, June 27, 2016 2:12 AM > *To:* firedrake > *Subject:* Re: [firedrake] Stochastic differential equations > > Hi! > > Yes, you can just use the numpy random number generator to write to the > basis coefficients for this. > > from numpy import random > > F = Function(V) > F.dat.data[:] = random.randn(N) > > etc. > > The main thing you have to be careful of is MPI parallelism. > > We have been using this to implement a stochastic QG equation for example. > > It might be nice to integrate random number generators into TSFL generated > code so that we can use them directly in UFL expressions - that would tidy > up a lot of our code (and help make it more parallel safe). > > cheers > --cjc > > On 26 June 2016 at 21:27, Angwenyi David <kipkoej(a)gmail.com> wrote: > >> Thanks, Francis. >> >> True, the DE is driven by white noise, that is, I assume it is gaussian >> with zero mean vector and identity covariance matrix premultiplied by the >> square root of the time-step. >> >> On 26 Jun 2016, at 22:16, Francis Poulin <fpoulin(a)uwaterloo.ca> wrote: >> >> Hello, >> >> I can't say whether there is anything specific to firedrake but I can say >> that in general it depends on what kind of noise you are considering in the >> SDE. >> >> If you are using coloured noise, then you treat your forcing term as you >> would any other term. However, given the fact that you are using W I >> imagine you are using White noise. >> >> In that case the popular thing to do is to use a splitting method. For >> the "regular" deterministic terms, you use whatever time stepping scheme >> you normally would use. However, for the white noise you probably want to >> use an Euler step. Basically, that means you have square root of dt in >> front of the random variable. The stochastic part is only first order >> accurate but it's much easier than using other techniques, but they do >> exist. >> >> Anyhow, that's my two cents worth. >> >> Francis >> >> ------------------ >> Francis Poulin >> Associate Professor >> Department of Applied Mathematics >> University of Waterloo >> >> email: fpoulin(a)uwaterloo.ca >> Web: https://uwaterloo.ca/poulin-research-group/ >> Telephone: +1 519 888 4567 x32637 >> >> ------------------------------ >> *From:* firedrake-bounces(a)imperial.ac.uk [ >> firedrake-bounces(a)imperial.ac.uk] on behalf of Angwenyi David [ >> kipkoej(a)gmail.com] >> *Sent:* Sunday, June 26, 2016 4:11 PM >> *To:* firedrake(a)imperial.ac.uk >> *Subject:* [firedrake] Stochastic differential equations >> >> Dear all, >> >> Is it possible to solve pdes driven by some noise using Firedrake? For >> example, see equation 3.2 in the attachment whose discretisation yields a >> system 3.3 and 3.4, where dW(t):=sqrt(dt)N(0,Identity). My other related >> concern is whether It is possible to implement, say, A Kalman filter with >> equation 3.2 as model equation on Firedrake. >> >> Thanks in advance. >> >> Angwenyi David >> <Screen Shot 2016-06-26 at 21.51.52.png> >> _______________________________________________ >> firedrake mailing list >> firedrake(a)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

1 0

Re: [firedrake] problems updating and running firedrake
by David Ham 27 Jun '16

27 Jun '16

Hi Francis, It looks like you might have hit the issue they describe here: http://askubuntu.com/questions/509283/python-no-module-named-datetime Basically datetime was apparently rolled into the python binary in a 2.7 point release. I don't know why this was changed but it makes the python binary in your virtualenv incompatible with your system python install. According to the thread on askubuntu, you have two options. a) Apparently you can re-initialise the virtualenv by running: virtualenv . in your top Firedrake directory. I'm guessing you should do this with the virtualenv deactivated. or b) A fresh Firedrake install should work because it will create a fresh virtualenv based on current Python binary. Regards, David On Mon, 27 Jun 2016 at 15:20 Francis Poulin <fpoulin(a)uwaterloo.ca> wrote: > Hello, > > I recently updated my ubuntu from 14.04 to 16.04. I tried running > firedrake and it failed. > > I also tried updating firedrake and it failed, see below. > > Do I need to install it from scratch or should it be possible to update it > with a new OS? > > Cheers, Francis > > > > Installing pip dependencies for COFFEE > /home/fpoulin/software/firedrake/bin/pip install -r COFFEE/requirements.txt > > Traceback (most recent call last): > File "/home/fpoulin/software/firedrake/bin/pip", line 7, in <module> > from pip import main > File > "/home/fpoulin/software/firedrake/local/lib/python2.7/site-packages/pip/__init__.py", > line 14, in <module> > from pip.utils import get_installed_distributions, get_prog > File > "/home/fpoulin/software/firedrake/local/lib/python2.7/site-packages/pip/utils/__init__.py", > line 27, in <module> > from pip._vendor import pkg_resources > File > "/home/fpoulin/software/firedrake/local/lib/python2.7/site-packages/pip/_vendor/pkg_resources/__init__.py", > line 36, in <module> > import plistlib > File "/usr/lib/python2.7/plistlib.py", line 62, in <module> > import datetime > ImportError: No module named datetime > Traceback (most recent call last): > File "../bin/firedrake-update", line 974, in <module> > pip_requirements(p) > File "../bin/firedrake-update", line 377, in pip_requirements > run_pip_install(["-r", "%s/requirements.txt" % package]) > File "../bin/firedrake-update", line 364, in run_pip_install > check_call(pipinstall + pipargs) > File "../bin/firedrake-update", line 211, in check_call > subprocess.check_call(arguments, env=env) > File "/usr/lib/python2.7/subprocess.py", line 541, in check_call > raise CalledProcessError(retcode, cmd) > subprocess.CalledProcessError: Command > '['/home/fpoulin/software/firedrake/bin/pip', 'install', '-r', > 'COFFEE/requirements.txt']' returned non-zero exit status 1 > > ------------------ > Francis Poulin > Associate Professor > Department of Applied Mathematics > University of Waterloo > > email: fpoulin(a)uwaterloo.ca > Web: https://uwaterloo.ca/poulin-research-group/ > Telephone: +1 519 888 4567 x32637 > >

2 1

problems updating and running firedrake
by Francis Poulin 27 Jun '16

27 Jun '16

Hello, I recently updated my ubuntu from 14.04 to 16.04. I tried running firedrake and it failed. I also tried updating firedrake and it failed, see below. Do I need to install it from scratch or should it be possible to update it with a new OS? Cheers, Francis Installing pip dependencies for COFFEE /home/fpoulin/software/firedrake/bin/pip install -r COFFEE/requirements.txt Traceback (most recent call last): File "/home/fpoulin/software/firedrake/bin/pip", line 7, in <module> from pip import main File "/home/fpoulin/software/firedrake/local/lib/python2.7/site-packages/pip/__init__.py", line 14, in <module> from pip.utils import get_installed_distributions, get_prog File "/home/fpoulin/software/firedrake/local/lib/python2.7/site-packages/pip/utils/__init__.py", line 27, in <module> from pip._vendor import pkg_resources File "/home/fpoulin/software/firedrake/local/lib/python2.7/site-packages/pip/_vendor/pkg_resources/__init__.py", line 36, in <module> import plistlib File "/usr/lib/python2.7/plistlib.py", line 62, in <module> import datetime ImportError: No module named datetime Traceback (most recent call last): File "../bin/firedrake-update", line 974, in <module> pip_requirements(p) File "../bin/firedrake-update", line 377, in pip_requirements run_pip_install(["-r", "%s/requirements.txt" % package]) File "../bin/firedrake-update", line 364, in run_pip_install check_call(pipinstall + pipargs) File "../bin/firedrake-update", line 211, in check_call subprocess.check_call(arguments, env=env) File "/usr/lib/python2.7/subprocess.py", line 541, in check_call raise CalledProcessError(retcode, cmd) subprocess.CalledProcessError: Command '['/home/fpoulin/software/firedrake/bin/pip', 'install', '-r', 'COFFEE/requirements.txt']' returned non-zero exit status 1 ------------------ Francis Poulin Associate Professor Department of Applied Mathematics University of Waterloo email: fpoulin(a)uwaterloo.ca Web: https://uwaterloo.ca/poulin-research-group/ Telephone: +1 519 888 4567 x32637

1 0

Re: [firedrake] Stochastic differential equations
by Colin Cotter 27 Jun '16

27 Jun '16

Hi! Yes, you can just use the numpy random number generator to write to the basis coefficients for this. from numpy import random F = Function(V) F.dat.data[:] = random.randn(N) etc. The main thing you have to be careful of is MPI parallelism. We have been using this to implement a stochastic QG equation for example. It might be nice to integrate random number generators into TSFL generated code so that we can use them directly in UFL expressions - that would tidy up a lot of our code (and help make it more parallel safe). cheers --cjc On 26 June 2016 at 21:27, Angwenyi David <kipkoej(a)gmail.com> wrote: > Thanks, Francis. > > True, the DE is driven by white noise, that is, I assume it is gaussian > with zero mean vector and identity covariance matrix premultiplied by the > square root of the time-step. > > On 26 Jun 2016, at 22:16, Francis Poulin <fpoulin(a)uwaterloo.ca> wrote: > > Hello, > > I can't say whether there is anything specific to firedrake but I can say > that in general it depends on what kind of noise you are considering in the > SDE. > > If you are using coloured noise, then you treat your forcing term as you > would any other term. However, given the fact that you are using W I > imagine you are using White noise. > > In that case the popular thing to do is to use a splitting method. For > the "regular" deterministic terms, you use whatever time stepping scheme > you normally would use. However, for the white noise you probably want to > use an Euler step. Basically, that means you have square root of dt in > front of the random variable. The stochastic part is only first order > accurate but it's much easier than using other techniques, but they do > exist. > > Anyhow, that's my two cents worth. > > Francis > > ------------------ > Francis Poulin > Associate Professor > Department of Applied Mathematics > University of Waterloo > > email: fpoulin(a)uwaterloo.ca > Web: https://uwaterloo.ca/poulin-research-group/ > Telephone: +1 519 888 4567 x32637 > > ------------------------------ > *From:* firedrake-bounces(a)imperial.ac.uk [firedrake-bounces(a)imperial.ac.uk] > on behalf of Angwenyi David [kipkoej(a)gmail.com] > *Sent:* Sunday, June 26, 2016 4:11 PM > *To:* firedrake(a)imperial.ac.uk > *Subject:* [firedrake] Stochastic differential equations > > Dear all, > > Is it possible to solve pdes driven by some noise using Firedrake? For > example, see equation 3.2 in the attachment whose discretisation yields a > system 3.3 and 3.4, where dW(t):=sqrt(dt)N(0,Identity). My other related > concern is whether It is possible to implement, say, A Kalman filter with > equation 3.2 as model equation on Firedrake. > > Thanks in advance. > > Angwenyi David > <Screen Shot 2016-06-26 at 21.51.52.png> > _______________________________________________ > firedrake mailing list > firedrake(a)imperial.ac.uk > https://mailman.ic.ac.uk/mailman/listinfo/firedrake > > > -- http://www.imperial.ac.uk/people/colin.cotter www.cambridge.org/9781107663916

3 3

Re: [firedrake] Stochastic differential equations
by Colin Cotter 27 Jun '16

27 Jun '16

Dear Angwenyi, For an ensemble filtering approach, please can you talk to my student Ali Gregory about how he did it? http://www.imperial.ac.uk/people/a.gregory14 If you are working with an ensemble, then one approach is to use a list of Firedrake fields for the different states. Then, the Kalman update needs to be implemented through a PyOP2 parloop. Ali can help you with this (he has done it for the ETPF). all the best --cjc On 27 June 2016 at 12:16, Angwenyi David <kipkoej(a)gmail.com> wrote: > The numpy random number generator works! Thanks Colin. > > Angwenyi > > On 27 Jun 2016, at 09:16, Angwenyi David <kipkoej(a)gmail.com> wrote: > > Thanks, Collins. > > Could you please, if it won’t matter to you, send me the code for the > stochastic quasi geostrophic equation you implemented. It shall help me a > great deal. > > How about implementing a Kalman filter with the system of equations 3.3 > and 3.4 as a model equations with given data? I presume it is possible; for > the system can be represented as psi^(n+1)=A^(n)*psi^(n)+noise. The > complication I see is having the spatial discretisation of phi fit into the > matrix A. Links to the articles with similar problems are appreciated. > > Thanks. > > Angwenyi > > On 27 Jun 2016, at 08:12, Colin Cotter <colin.cotter(a)imperial.ac.uk> > wrote: > > Hi! > > Yes, you can just use the numpy random number generator to write to the > basis coefficients for this. > > from numpy import random > > F = Function(V) > F.dat.data[:] = random.randn(N) > > etc. > > The main thing you have to be careful of is MPI parallelism. > > We have been using this to implement a stochastic QG equation for example. > > It might be nice to integrate random number generators into TSFL generated > code so that we can use them directly in UFL expressions - that would tidy > up a lot of our code (and help make it more parallel safe). > > cheers > --cjc > > On 26 June 2016 at 21:27, Angwenyi David <kipkoej(a)gmail.com> wrote: > >> Thanks, Francis. >> >> True, the DE is driven by white noise, that is, I assume it is gaussian >> with zero mean vector and identity covariance matrix premultiplied by the >> square root of the time-step. >> >> On 26 Jun 2016, at 22:16, Francis Poulin <fpoulin(a)uwaterloo.ca> wrote: >> >> Hello, >> >> I can't say whether there is anything specific to firedrake but I can say >> that in general it depends on what kind of noise you are considering in the >> SDE. >> >> If you are using coloured noise, then you treat your forcing term as you >> would any other term. However, given the fact that you are using W I >> imagine you are using White noise. >> >> In that case the popular thing to do is to use a splitting method. For >> the "regular" deterministic terms, you use whatever time stepping scheme >> you normally would use. However, for the white noise you probably want to >> use an Euler step. Basically, that means you have square root of dt in >> front of the random variable. The stochastic part is only first order >> accurate but it's much easier than using other techniques, but they do >> exist. >> >> Anyhow, that's my two cents worth. >> >> Francis >> >> ------------------ >> Francis Poulin >> Associate Professor >> Department of Applied Mathematics >> University of Waterloo >> >> email: fpoulin(a)uwaterloo.ca >> Web: https://uwaterloo.ca/poulin-research-group/ >> Telephone: +1 519 888 4567 x32637 >> >> ------------------------------ >> *From:* firedrake-bounces(a)imperial.ac.uk [ >> firedrake-bounces(a)imperial.ac.uk] on behalf of Angwenyi David [ >> kipkoej(a)gmail.com] >> *Sent:* Sunday, June 26, 2016 4:11 PM >> *To:* firedrake(a)imperial.ac.uk >> *Subject:* [firedrake] Stochastic differential equations >> >> Dear all, >> >> Is it possible to solve pdes driven by some noise using Firedrake? For >> example, see equation 3.2 in the attachment whose discretisation yields a >> system 3.3 and 3.4, where dW(t):=sqrt(dt)N(0,Identity). My other related >> concern is whether It is possible to implement, say, A Kalman filter with >> equation 3.2 as model equation on Firedrake. >> >> Thanks in advance. >> >> Angwenyi David >> <Screen Shot 2016-06-26 at 21.51.52.png> >> _______________________________________________ >> firedrake mailing list >> firedrake(a)imperial.ac.uk >> https://mailman.ic.ac.uk/mailman/listinfo/firedrake >> >> >> > > > -- > http://www.imperial.ac.uk/people/colin.cotter > > www.cambridge.org/9781107663916 > > > _______________________________________________ > firedrake mailing list > firedrake(a)imperial.ac.uk > https://mailman.ic.ac.uk/mailman/listinfo/firedrake > > > > -- http://www.imperial.ac.uk/people/colin.cotter www.cambridge.org/9781107663916

2 1

Re: [firedrake] upwind example remark; extension to CN in time
by Andrew McRae 27 Jun '16

27 Jun '16

(For simplicity, I'll assume an explicit timestepping method) You have a field from your previous timestep; let's call this U_ You represent your boundary data using some other field; let's call this U0. (Note: Firedrake can't yet represent fields that "only exist on the boundary", so U0 will just be a normal global function, but we will only make use of its boundary values, so it can be nonsense or zero in the interior) On the LHS of your equation, you just have a mass matrix (assume you are solving for the new field U, or for the increment dU, it doesn't matter) On the RHS, you have some cell integral after integrating by parts, some integral over interior facets that uses U_('+'), U_('-'), n('+') and n('-')***, and some integral over exterior facets that uses U_, U0, and n. If your interior facet flux integral was dot(jump(phi), un('+')*U('+') - un( '-')*U('-'))*dS (where un is u.n if this is positive, else 0; u is the velocity field), your exterior facet flux integral could be dot(phi, un*U_ + (dot(u, n) - un)*U0)*ds [Gevaar: may contain sign errors] ***: Usually n('-') is -n('+'), but since we separate the mesh geometry from the mesh topology, this doesn't have to be true <http://www.firedrakeproject.org/extruded-meshes.html#hedgehog-extrusion>. On 27 June 2016 at 11:40, Onno Bokhove <O.Bokhove(a)leeds.ac.uk> wrote: > So how would you implement a standard numerical flux at the > > boundary that uses a mix of information from the outside and given > > -which I would call "+", as well as the interior or "-" information? > > In other words it needs a U+ and a U-? > > > Using the standard and same numerical flux used everywhere else > > in the interior (i.e. needing a U+ and U- as input as well) except that > one piece of information is supplied? Basically a ghost value? Is > that not allowed? > > > > ------------------------------ > *From:* firedrake-bounces(a)imperial.ac.uk <firedrake-bounces(a)imperial.ac.uk> > on behalf of Andrew McRae <A.T.T.McRae(a)bath.ac.uk> > *Sent:* Monday, June 27, 2016 11:33 AM > > *To:* firedrake(a)imperial.ac.uk > *Subject:* Re: [firedrake] upwind example remark; extension to CN in time > > Again, we do not use + and - on exterior facets. > > On 27 June 2016 at 11:29, Onno Bokhove <O.Bokhove(a)leeds.ac.uk> wrote: > >> Ah, that is: in our usual set-up at any face there is an orientation >> chosen with one of the normals being called outside. >> >> That is arbitrary but once fixed it is known. >> >> In the interior that does not matter but at real boundaries >> >> the convention is always that the normal assigned is the >> >> normal pointing outside the domain. >> >> That is why at real boundaries one knows which sign convention >> >> concerns outside information. >> >> Otherwise how do you know whether to assign an outside the >> >> domain value to + or -? You would not know. >> >> >> In my convention with - inside, the + side is then the given >> >> information from the outside. >> >> >> ------------------------------ >> *From:* firedrake-bounces(a)imperial.ac.uk < >> firedrake-bounces(a)imperial.ac.uk> on behalf of Colin Cotter < >> colin.cotter(a)imperial.ac.uk> >> *Sent:* Monday, June 27, 2016 7:41 AM >> *To:* firedrake >> *Subject:* Re: [firedrake] upwind example remark; extension to CN in time >> >> Hi Onno, >> But + and - are not on the inside or outside, they are just two >> arbitrary sides. >> >> I'll look into why the code is crashing today hopefully, I have a phone >> conference about admissions until 4pm. >> >> all the best >> >> --cjc >> >> On 27 June 2016 at 04:34, Onno Bokhove <O.Bokhove(a)leeds.ac.uk> wrote: >> >>> Hi Colin, >>> >>> >>> I know that of course as that is how I loop in my own c or matlab >>> >>> programs! One loop over elements and one loop over the faces. >>> >>> >>> I am just saying state in the documentation example >>> >>> that plus is on the inside when the normal is pointing outward, >>> >>> and minus on the outside. Just add one sentence. >>> >>> I usually tend to do that the other way around as >>> >>> it is the limit with -epsilon and epsilon positive that is on the inside >>> (so minus not plus is my usual convention; it does not matter but define it >>> also >>> >>> in the documentation example for clarity). >>> >>> >>> * The one remaining issue in my code seems that FD does not like me to >>> use D_(+) or D_(-) while D(+) and D(-) are allowed, with the one >>> >>> difference being that D is a TrailFunction >>> >>> and D_ (previous time step) a Function. >>> >>> I guess I don't understand that difference yet. >>> >>> >>> Thanks for the UFL link. That will help. >>> >>> >>> Best wishes, >>> >>> >>> >>> >>> >>> >>> ------------------------------ >>> *From:* firedrake-bounces(a)imperial.ac.uk < >>> firedrake-bounces(a)imperial.ac.uk> on behalf of Colin Cotter < >>> colin.cotter(a)imperial.ac.uk> >>> *Sent:* Sunday, June 26, 2016 3:11 PM >>> >>> *To:* firedrake >>> *Subject:* Re: [firedrake] upwind example remark; extension to CN in >>> time >>> >>> Hi Onno, >>> This is documented in the UFL manual (the plusses, minusses etc.). >>> >>> fenicsproject.org/pub/documents/*ufl*/*ufl*-user-*manual*/*ufl*-user- >>> *manual*.pdf >>> >>> The reason for this notation is that UFL is written for a global domain, >>> in particular for facets it is the union of all interior facets, and when >>> surface integrals are done, the facets are visited one by one instead of >>> the cells. This means that there needs to be an (arbitrary) labelling of + >>> and - sides of each facet. The notation was originally invented to write >>> down global inner products to analyse DG methods for elliptic problems, >>> where it is completely standard. This notation was set out in >>> Arnold, Douglas N., et al. "Unified analysis of discontinuous Galerkin >>> methods for elliptic problems." *SIAM journal on numerical analysis* 39.5 >>> (2002): 1749-1779. >>> >>> all the best >>> --cjc >>> >>> On 26 June 2016 at 14:08, Onno Bokhove <O.Bokhove(a)leeds.ac.uk> wrote: >>> >>>> Sure I knew that. D+ is on the inside of a cell with n pointing outward >>>> >>>> unless I am blinded by being on the other side of the world. >>>> >>>> Bit unusual convention in the DG world. Usually minus is on the inside. >>>> >>>> But could it as such be defined in the write-up/documentation please? >>>> >>>> In any DG article the convention whatever it is gets defined. >>>> >>>> >>>> ------------------------------ >>>> *From:* firedrake-bounces(a)imperial.ac.uk < >>>> firedrake-bounces(a)imperial.ac.uk> on behalf of Andrew McRae < >>>> A.T.T.McRae(a)bath.ac.uk> >>>> *Sent:* Sunday, June 26, 2016 12:40 PM >>>> *To:* firedrake(a)imperial.ac.uk >>>> >>>> *Subject:* Re: [firedrake] upwind example remark; extension to CN in >>>> time >>>> >>>> D+ and D- refer to the values either side of an interior facet between >>>> two cells. >>>> >>>> On the boundary (exterior facets), there is no such ambiguity, so just >>>> D is used. >>>> >>>> In the steady-state example on the website, this D enters the >>>> left-hand-side integral over the outflow part of the boundary, while the >>>> prescribed function D0 is used in the integral over the inflow part. >>>> >>>> "In general", you can write a term that contains both D and D0 and >>>> 'picks the right one' depending on the sign of u.n (this is how the >>>> interior facet integral works). Note that n always points *outwards* on >>>> the exterior boundary facets. >>>> >>>> On 26 June 2016 at 12:24, Onno Bokhove <O.Bokhove(a)leeds.ac.uk> wrote: >>>> >>>>> The issue is that D_('+') (or D_('-')) is not defined but I don't >>>>> understand why not, or how to fix it? >>>>> >>>>> >>>>> ------------------------------ >>>>> *From:* firedrake-bounces(a)imperial.ac.uk < >>>>> firedrake-bounces(a)imperial.ac.uk> on behalf of Onno Bokhove < >>>>> O.Bokhove(a)leeds.ac.uk> >>>>> *Sent:* Sunday, June 26, 2016 12:21 PM >>>>> >>>>> *To:* firedrake >>>>> *Subject:* Re: [firedrake] upwind example remark; extension to CN in >>>>> time >>>>> >>>>> >>>>> Yes that was my point indeed. >>>>> >>>>> To build any more advanced code upon this one, one >>>>> >>>>> would need to use an alternate approach entirely >>>>> >>>>> through the flux. >>>>> >>>>> That is why it may be useful to extend the >>>>> >>>>> example with two alternatives: the simple one shown >>>>> >>>>> and one flux-based only. >>>>> >>>>> >>>>> By the way, in the attached I have L2=0 and thet=0; >>>>> >>>>> if I uncomment the actual L2 there is a complaint, >>>>> >>>>> even though they=0. >>>>> >>>>> Obviously I need this term for the general theta-scheme, >>>>> >>>>> including thet=0.5 (Crank-Nicolson). >>>>> >>>>> I am doing the above simple extensions to get used to FD! >>>>> Small steps ... >>>>> >>>>> The example is ties such that value 1 advects into through the >>>>> entire domain in the time set. >>>>> ------------------------------ >>>>> *From:* firedrake-bounces(a)imperial.ac.uk < >>>>> firedrake-bounces(a)imperial.ac.uk> on behalf of Colin Cotter < >>>>> colin.cotter(a)imperial.ac.uk> >>>>> *Sent:* Sunday, June 26, 2016 11:26 AM >>>>> *To:* firedrake >>>>> *Subject:* Re: [firedrake] upwind example remark; extension to CN in >>>>> time >>>>> >>>>> Hi Onno, >>>>> In answer to your first question. Things are simplified in the >>>>> advection demo because we know which edges are upwind, which are downwind, >>>>> and which have no flow through them. If you wanted to generalise this, then >>>>> we would need a surface integral over the entire external boundary with a >>>>> switch that either gave the internal value of the tracer or the inflow >>>>> boundary condition, depending on the direction of u.n, just like we did for >>>>> the interior facets. >>>>> Instead, in this example, on the outflow boundary we just write an >>>>> integral involving the interior value, and on the inflow boundary we write >>>>> an integral setting the inflow values. >>>>> >>>>> all the best >>>>> -cjc >>>>> >>>>> On 25 June 2016 at 16:05, Onno Bokhove <O.Bokhove(a)leeds.ac.uk> wrote: >>>>> >>>>>> Hi FDs, >>>>>> >>>>>> >>>>>> (i) The upwind example is odd in that the real boundaries are treated >>>>>> >>>>>> in a special way rather than setting the outside value in >>>>>> >>>>>> the same numerical flux as anywhere else; >>>>>> >>>>>> my question is how this can be fixed? >>>>>> >>>>>> So no special terms at the in- and outflow, just specify >>>>>> >>>>>> the outside velocity used. This is imo the proper DG way not was is >>>>>> done >>>>>> >>>>>> now. >>>>>> >>>>>> >>>>>> (ii) Attached my attempt to do a time dependent advection >>>>>> >>>>>> with Crank Nicolson, which works but does not yet >>>>>> >>>>>> give the right answer, in part due to (i). >>>>>> >>>>>> >>>>>> [I don't understand why I have to make >>>>>> >>>>>> my theta variable that a constant but other wise it did not work.] >>>>>> >>>>>> Any suggestions, corrections, also on other silly things done by me? >>>>>> >>>>>> Next step is to extend to system with 3 equations; >>>>>> hence the few extra but unused variables. >>>>>> >>>>>> Thanks, from midnight Seoul in old style Korean guesthouse. >>>>>> >>>>>> Onno >>>>>> ------------------------------ >>>>>> >>>>> >>>>> >>>>> >>>>> -- >>>>> http://www.imperial.ac.uk/people/colin.cotter >>>>> >>>>> www.cambridge.org/9781107663916 >>>>> >>>>> >>>>> >>>> >>> >>> >>> -- >>> http://www.imperial.ac.uk/people/colin.cotter >>> >>> www.cambridge.org/9781107663916 >>> >>> >>> >> >> >> -- >> http://www.imperial.ac.uk/people/colin.cotter >> >> www.cambridge.org/9781107663916 >> >> >> >

1 0

Re: [firedrake] upwind example remark; extension to CN in time
by David Ham 27 Jun '16

27 Jun '16

On an exterior facet, FacetNormal points out of the domain. This is sufficient. For example, suppose I have a velocity field u, then: u_out = 0.5* (dot(u, n) + abs(dot(u, n))) is equal to the outward flux when the velocity goes out of the domain, and 0 when the velocity is coming in. This is sufficient to select between the field value inside the domain, and a prescribed boundary value. More complex flux decisions can be created in a similar manner. Regards, David On Mon, 27 Jun 2016 at 11:43 Onno Bokhove <O.Bokhove(a)leeds.ac.uk> wrote: > So how would you implement a standard numerical flux at the > > boundary that uses a mix of information from the outside and given > > -which I would call "+", as well as the interior or "-" information? > > In other words it needs a U+ and a U-? > > > Using the standard and same numerical flux used everywhere else > > in the interior (i.e. needing a U+ and U- as input as well) except that > one piece of information is supplied? Basically a ghost value? Is > that not allowed? > > > > ------------------------------ > *From:* firedrake-bounces(a)imperial.ac.uk <firedrake-bounces(a)imperial.ac.uk> > on behalf of Andrew McRae <A.T.T.McRae(a)bath.ac.uk> > *Sent:* Monday, June 27, 2016 11:33 AM > > *To:* firedrake(a)imperial.ac.uk > *Subject:* Re: [firedrake] upwind example remark; extension to CN in time > Again, we do not use + and - on exterior facets. > > On 27 June 2016 at 11:29, Onno Bokhove <O.Bokhove(a)leeds.ac.uk> wrote: > >> Ah, that is: in our usual set-up at any face there is an orientation >> chosen with one of the normals being called outside. >> >> That is arbitrary but once fixed it is known. >> >> In the interior that does not matter but at real boundaries >> >> the convention is always that the normal assigned is the >> >> normal pointing outside the domain. >> >> That is why at real boundaries one knows which sign convention >> >> concerns outside information. >> >> Otherwise how do you know whether to assign an outside the >> >> domain value to + or -? You would not know. >> >> >> In my convention with - inside, the + side is then the given >> >> information from the outside. >> >> >> ------------------------------ >> *From:* firedrake-bounces(a)imperial.ac.uk < >> firedrake-bounces(a)imperial.ac.uk> on behalf of Colin Cotter < >> colin.cotter(a)imperial.ac.uk> >> *Sent:* Monday, June 27, 2016 7:41 AM >> *To:* firedrake >> *Subject:* Re: [firedrake] upwind example remark; extension to CN in time >> >> Hi Onno, >> But + and - are not on the inside or outside, they are just two >> arbitrary sides. >> >> I'll look into why the code is crashing today hopefully, I have a phone >> conference about admissions until 4pm. >> >> all the best >> >> --cjc >> >> On 27 June 2016 at 04:34, Onno Bokhove <O.Bokhove(a)leeds.ac.uk> wrote: >> >>> Hi Colin, >>> >>> >>> I know that of course as that is how I loop in my own c or matlab >>> >>> programs! One loop over elements and one loop over the faces. >>> >>> >>> I am just saying state in the documentation example >>> >>> that plus is on the inside when the normal is pointing outward, >>> >>> and minus on the outside. Just add one sentence. >>> >>> I usually tend to do that the other way around as >>> >>> it is the limit with -epsilon and epsilon positive that is on the inside >>> (so minus not plus is my usual convention; it does not matter but define it >>> also >>> >>> in the documentation example for clarity). >>> >>> >>> * The one remaining issue in my code seems that FD does not like me to >>> use D_(+) or D_(-) while D(+) and D(-) are allowed, with the one >>> >>> difference being that D is a TrailFunction >>> >>> and D_ (previous time step) a Function. >>> >>> I guess I don't understand that difference yet. >>> >>> >>> Thanks for the UFL link. That will help. >>> >>> >>> Best wishes, >>> >>> >>> >>> >>> >>> >>> ------------------------------ >>> *From:* firedrake-bounces(a)imperial.ac.uk < >>> firedrake-bounces(a)imperial.ac.uk> on behalf of Colin Cotter < >>> colin.cotter(a)imperial.ac.uk> >>> *Sent:* Sunday, June 26, 2016 3:11 PM >>> >>> *To:* firedrake >>> *Subject:* Re: [firedrake] upwind example remark; extension to CN in >>> time >>> >>> Hi Onno, >>> This is documented in the UFL manual (the plusses, minusses etc.). >>> >>> fenicsproject.org/pub/documents/*ufl*/*ufl*-user-*manual*/*ufl*-user- >>> *manual*.pdf >>> >>> The reason for this notation is that UFL is written for a global domain, >>> in particular for facets it is the union of all interior facets, and when >>> surface integrals are done, the facets are visited one by one instead of >>> the cells. This means that there needs to be an (arbitrary) labelling of + >>> and - sides of each facet. The notation was originally invented to write >>> down global inner products to analyse DG methods for elliptic problems, >>> where it is completely standard. This notation was set out in >>> Arnold, Douglas N., et al. "Unified analysis of discontinuous Galerkin >>> methods for elliptic problems." *SIAM journal on numerical analysis* 39.5 >>> (2002): 1749-1779. >>> >>> all the best >>> --cjc >>> >>> On 26 June 2016 at 14:08, Onno Bokhove <O.Bokhove(a)leeds.ac.uk> wrote: >>> >>>> Sure I knew that. D+ is on the inside of a cell with n pointing outward >>>> >>>> unless I am blinded by being on the other side of the world. >>>> >>>> Bit unusual convention in the DG world. Usually minus is on the inside. >>>> >>>> But could it as such be defined in the write-up/documentation please? >>>> >>>> In any DG article the convention whatever it is gets defined. >>>> >>>> >>>> ------------------------------ >>>> *From:* firedrake-bounces(a)imperial.ac.uk < >>>> firedrake-bounces(a)imperial.ac.uk> on behalf of Andrew McRae < >>>> A.T.T.McRae(a)bath.ac.uk> >>>> *Sent:* Sunday, June 26, 2016 12:40 PM >>>> *To:* firedrake(a)imperial.ac.uk >>>> >>>> *Subject:* Re: [firedrake] upwind example remark; extension to CN in >>>> time >>>> >>>> D+ and D- refer to the values either side of an interior facet between >>>> two cells. >>>> >>>> On the boundary (exterior facets), there is no such ambiguity, so just >>>> D is used. >>>> >>>> In the steady-state example on the website, this D enters the >>>> left-hand-side integral over the outflow part of the boundary, while the >>>> prescribed function D0 is used in the integral over the inflow part. >>>> >>>> "In general", you can write a term that contains both D and D0 and >>>> 'picks the right one' depending on the sign of u.n (this is how the >>>> interior facet integral works). Note that n always points *outwards* on >>>> the exterior boundary facets. >>>> >>>> On 26 June 2016 at 12:24, Onno Bokhove <O.Bokhove(a)leeds.ac.uk> wrote: >>>> >>>>> The issue is that D_('+') (or D_('-')) is not defined but I don't >>>>> understand why not, or how to fix it? >>>>> >>>>> >>>>> ------------------------------ >>>>> *From:* firedrake-bounces(a)imperial.ac.uk < >>>>> firedrake-bounces(a)imperial.ac.uk> on behalf of Onno Bokhove < >>>>> O.Bokhove(a)leeds.ac.uk> >>>>> *Sent:* Sunday, June 26, 2016 12:21 PM >>>>> >>>>> *To:* firedrake >>>>> *Subject:* Re: [firedrake] upwind example remark; extension to CN in >>>>> time >>>>> >>>>> >>>>> Yes that was my point indeed. >>>>> >>>>> To build any more advanced code upon this one, one >>>>> >>>>> would need to use an alternate approach entirely >>>>> >>>>> through the flux. >>>>> >>>>> That is why it may be useful to extend the >>>>> >>>>> example with two alternatives: the simple one shown >>>>> >>>>> and one flux-based only. >>>>> >>>>> >>>>> By the way, in the attached I have L2=0 and thet=0; >>>>> >>>>> if I uncomment the actual L2 there is a complaint, >>>>> >>>>> even though they=0. >>>>> >>>>> Obviously I need this term for the general theta-scheme, >>>>> >>>>> including thet=0.5 (Crank-Nicolson). >>>>> >>>>> I am doing the above simple extensions to get used to FD! >>>>> Small steps ... >>>>> >>>>> The example is ties such that value 1 advects into through the >>>>> entire domain in the time set. >>>>> ------------------------------ >>>>> *From:* firedrake-bounces(a)imperial.ac.uk < >>>>> firedrake-bounces(a)imperial.ac.uk> on behalf of Colin Cotter < >>>>> colin.cotter(a)imperial.ac.uk> >>>>> *Sent:* Sunday, June 26, 2016 11:26 AM >>>>> *To:* firedrake >>>>> *Subject:* Re: [firedrake] upwind example remark; extension to CN in >>>>> time >>>>> >>>>> Hi Onno, >>>>> In answer to your first question. Things are simplified in the >>>>> advection demo because we know which edges are upwind, which are downwind, >>>>> and which have no flow through them. If you wanted to generalise this, then >>>>> we would need a surface integral over the entire external boundary with a >>>>> switch that either gave the internal value of the tracer or the inflow >>>>> boundary condition, depending on the direction of u.n, just like we did for >>>>> the interior facets. >>>>> Instead, in this example, on the outflow boundary we just write an >>>>> integral involving the interior value, and on the inflow boundary we write >>>>> an integral setting the inflow values. >>>>> >>>>> all the best >>>>> -cjc >>>>> >>>>> On 25 June 2016 at 16:05, Onno Bokhove <O.Bokhove(a)leeds.ac.uk> wrote: >>>>> >>>>>> Hi FDs, >>>>>> >>>>>> >>>>>> (i) The upwind example is odd in that the real boundaries are treated >>>>>> >>>>>> in a special way rather than setting the outside value in >>>>>> >>>>>> the same numerical flux as anywhere else; >>>>>> >>>>>> my question is how this can be fixed? >>>>>> >>>>>> So no special terms at the in- and outflow, just specify >>>>>> >>>>>> the outside velocity used. This is imo the proper DG way not was is >>>>>> done >>>>>> >>>>>> now. >>>>>> >>>>>> >>>>>> (ii) Attached my attempt to do a time dependent advection >>>>>> >>>>>> with Crank Nicolson, which works but does not yet >>>>>> >>>>>> give the right answer, in part due to (i). >>>>>> >>>>>> >>>>>> [I don't understand why I have to make >>>>>> >>>>>> my theta variable that a constant but other wise it did not work.] >>>>>> >>>>>> Any suggestions, corrections, also on other silly things done by me? >>>>>> >>>>>> Next step is to extend to system with 3 equations; >>>>>> hence the few extra but unused variables. >>>>>> >>>>>> Thanks, from midnight Seoul in old style Korean guesthouse. >>>>>> >>>>>> Onno >>>>>> ------------------------------ >>>>>> >>>>> >>>>> >>>>> >>>>> -- >>>>> http://www.imperial.ac.uk/people/colin.cotter >>>>> >>>>> www.cambridge.org/9781107663916 >>>>> >>>>> >>>>> >>>> >>> >>> >>> -- >>> http://www.imperial.ac.uk/people/colin.cotter >>> >>> www.cambridge.org/9781107663916 >>> >>> >>> >> >> >> -- >> http://www.imperial.ac.uk/people/colin.cotter >> >> www.cambridge.org/9781107663916 >> >> >> >

1 0

Re: [firedrake] upwind example remark; extension to CN in time
by Andrew McRae 27 Jun '16

27 Jun '16

Again, we do not use + and - on exterior facets. On 27 June 2016 at 11:29, Onno Bokhove <O.Bokhove(a)leeds.ac.uk> wrote: > Ah, that is: in our usual set-up at any face there is an orientation > chosen with one of the normals being called outside. > > That is arbitrary but once fixed it is known. > > In the interior that does not matter but at real boundaries > > the convention is always that the normal assigned is the > > normal pointing outside the domain. > > That is why at real boundaries one knows which sign convention > > concerns outside information. > > Otherwise how do you know whether to assign an outside the > > domain value to + or -? You would not know. > > > In my convention with - inside, the + side is then the given > > information from the outside. > > > ------------------------------ > *From:* firedrake-bounces(a)imperial.ac.uk <firedrake-bounces(a)imperial.ac.uk> > on behalf of Colin Cotter <colin.cotter(a)imperial.ac.uk> > *Sent:* Monday, June 27, 2016 7:41 AM > *To:* firedrake > *Subject:* Re: [firedrake] upwind example remark; extension to CN in time > > Hi Onno, > But + and - are not on the inside or outside, they are just two > arbitrary sides. > > I'll look into why the code is crashing today hopefully, I have a phone > conference about admissions until 4pm. > > all the best > > --cjc > > On 27 June 2016 at 04:34, Onno Bokhove <O.Bokhove(a)leeds.ac.uk> wrote: > >> Hi Colin, >> >> >> I know that of course as that is how I loop in my own c or matlab >> >> programs! One loop over elements and one loop over the faces. >> >> >> I am just saying state in the documentation example >> >> that plus is on the inside when the normal is pointing outward, >> >> and minus on the outside. Just add one sentence. >> >> I usually tend to do that the other way around as >> >> it is the limit with -epsilon and epsilon positive that is on the inside >> (so minus not plus is my usual convention; it does not matter but define it >> also >> >> in the documentation example for clarity). >> >> >> * The one remaining issue in my code seems that FD does not like me to >> use D_(+) or D_(-) while D(+) and D(-) are allowed, with the one >> >> difference being that D is a TrailFunction >> >> and D_ (previous time step) a Function. >> >> I guess I don't understand that difference yet. >> >> >> Thanks for the UFL link. That will help. >> >> >> Best wishes, >> >> >> >> >> >> >> ------------------------------ >> *From:* firedrake-bounces(a)imperial.ac.uk < >> firedrake-bounces(a)imperial.ac.uk> on behalf of Colin Cotter < >> colin.cotter(a)imperial.ac.uk> >> *Sent:* Sunday, June 26, 2016 3:11 PM >> >> *To:* firedrake >> *Subject:* Re: [firedrake] upwind example remark; extension to CN in time >> >> Hi Onno, >> This is documented in the UFL manual (the plusses, minusses etc.). >> >> fenicsproject.org/pub/documents/*ufl*/*ufl*-user-*manual*/*ufl*-user- >> *manual*.pdf >> >> The reason for this notation is that UFL is written for a global domain, >> in particular for facets it is the union of all interior facets, and when >> surface integrals are done, the facets are visited one by one instead of >> the cells. This means that there needs to be an (arbitrary) labelling of + >> and - sides of each facet. The notation was originally invented to write >> down global inner products to analyse DG methods for elliptic problems, >> where it is completely standard. This notation was set out in >> Arnold, Douglas N., et al. "Unified analysis of discontinuous Galerkin >> methods for elliptic problems." *SIAM journal on numerical analysis* 39.5 >> (2002): 1749-1779. >> >> all the best >> --cjc >> >> On 26 June 2016 at 14:08, Onno Bokhove <O.Bokhove(a)leeds.ac.uk> wrote: >> >>> Sure I knew that. D+ is on the inside of a cell with n pointing outward >>> >>> unless I am blinded by being on the other side of the world. >>> >>> Bit unusual convention in the DG world. Usually minus is on the inside. >>> >>> But could it as such be defined in the write-up/documentation please? >>> >>> In any DG article the convention whatever it is gets defined. >>> >>> >>> ------------------------------ >>> *From:* firedrake-bounces(a)imperial.ac.uk < >>> firedrake-bounces(a)imperial.ac.uk> on behalf of Andrew McRae < >>> A.T.T.McRae(a)bath.ac.uk> >>> *Sent:* Sunday, June 26, 2016 12:40 PM >>> *To:* firedrake(a)imperial.ac.uk >>> >>> *Subject:* Re: [firedrake] upwind example remark; extension to CN in >>> time >>> >>> D+ and D- refer to the values either side of an interior facet between >>> two cells. >>> >>> On the boundary (exterior facets), there is no such ambiguity, so just D >>> is used. >>> >>> In the steady-state example on the website, this D enters the >>> left-hand-side integral over the outflow part of the boundary, while the >>> prescribed function D0 is used in the integral over the inflow part. >>> >>> "In general", you can write a term that contains both D and D0 and >>> 'picks the right one' depending on the sign of u.n (this is how the >>> interior facet integral works). Note that n always points *outwards* on >>> the exterior boundary facets. >>> >>> On 26 June 2016 at 12:24, Onno Bokhove <O.Bokhove(a)leeds.ac.uk> wrote: >>> >>>> The issue is that D_('+') (or D_('-')) is not defined but I don't >>>> understand why not, or how to fix it? >>>> >>>> >>>> ------------------------------ >>>> *From:* firedrake-bounces(a)imperial.ac.uk < >>>> firedrake-bounces(a)imperial.ac.uk> on behalf of Onno Bokhove < >>>> O.Bokhove(a)leeds.ac.uk> >>>> *Sent:* Sunday, June 26, 2016 12:21 PM >>>> >>>> *To:* firedrake >>>> *Subject:* Re: [firedrake] upwind example remark; extension to CN in >>>> time >>>> >>>> >>>> Yes that was my point indeed. >>>> >>>> To build any more advanced code upon this one, one >>>> >>>> would need to use an alternate approach entirely >>>> >>>> through the flux. >>>> >>>> That is why it may be useful to extend the >>>> >>>> example with two alternatives: the simple one shown >>>> >>>> and one flux-based only. >>>> >>>> >>>> By the way, in the attached I have L2=0 and thet=0; >>>> >>>> if I uncomment the actual L2 there is a complaint, >>>> >>>> even though they=0. >>>> >>>> Obviously I need this term for the general theta-scheme, >>>> >>>> including thet=0.5 (Crank-Nicolson). >>>> >>>> I am doing the above simple extensions to get used to FD! >>>> Small steps ... >>>> >>>> The example is ties such that value 1 advects into through the >>>> entire domain in the time set. >>>> ------------------------------ >>>> *From:* firedrake-bounces(a)imperial.ac.uk < >>>> firedrake-bounces(a)imperial.ac.uk> on behalf of Colin Cotter < >>>> colin.cotter(a)imperial.ac.uk> >>>> *Sent:* Sunday, June 26, 2016 11:26 AM >>>> *To:* firedrake >>>> *Subject:* Re: [firedrake] upwind example remark; extension to CN in >>>> time >>>> >>>> Hi Onno, >>>> In answer to your first question. Things are simplified in the >>>> advection demo because we know which edges are upwind, which are downwind, >>>> and which have no flow through them. If you wanted to generalise this, then >>>> we would need a surface integral over the entire external boundary with a >>>> switch that either gave the internal value of the tracer or the inflow >>>> boundary condition, depending on the direction of u.n, just like we did for >>>> the interior facets. >>>> Instead, in this example, on the outflow boundary we just write an >>>> integral involving the interior value, and on the inflow boundary we write >>>> an integral setting the inflow values. >>>> >>>> all the best >>>> -cjc >>>> >>>> On 25 June 2016 at 16:05, Onno Bokhove <O.Bokhove(a)leeds.ac.uk> wrote: >>>> >>>>> Hi FDs, >>>>> >>>>> >>>>> (i) The upwind example is odd in that the real boundaries are treated >>>>> >>>>> in a special way rather than setting the outside value in >>>>> >>>>> the same numerical flux as anywhere else; >>>>> >>>>> my question is how this can be fixed? >>>>> >>>>> So no special terms at the in- and outflow, just specify >>>>> >>>>> the outside velocity used. This is imo the proper DG way not was is >>>>> done >>>>> >>>>> now. >>>>> >>>>> >>>>> (ii) Attached my attempt to do a time dependent advection >>>>> >>>>> with Crank Nicolson, which works but does not yet >>>>> >>>>> give the right answer, in part due to (i). >>>>> >>>>> >>>>> [I don't understand why I have to make >>>>> >>>>> my theta variable that a constant but other wise it did not work.] >>>>> >>>>> Any suggestions, corrections, also on other silly things done by me? >>>>> >>>>> Next step is to extend to system with 3 equations; >>>>> hence the few extra but unused variables. >>>>> >>>>> Thanks, from midnight Seoul in old style Korean guesthouse. >>>>> >>>>> Onno >>>>> ------------------------------ >>>>> >>>> >>>> >>>> >>>> -- >>>> http://www.imperial.ac.uk/people/colin.cotter >>>> >>>> www.cambridge.org/9781107663916 >>>> >>>> >>>> >>> >> >> >> -- >> http://www.imperial.ac.uk/people/colin.cotter >> >> www.cambridge.org/9781107663916 >> >> >> > > > -- > http://www.imperial.ac.uk/people/colin.cotter > > www.cambridge.org/9781107663916 > > >

2 1

Re: [firedrake] upwind example remark; extension to CN in time
by Colin Cotter 27 Jun '16

27 Jun '16

Hi Onno, But + and - are not on the inside or outside, they are just two arbitrary sides. I'll look into why the code is crashing today hopefully, I have a phone conference about admissions until 4pm. all the best --cjc On 27 June 2016 at 04:34, Onno Bokhove <O.Bokhove(a)leeds.ac.uk> wrote: > Hi Colin, > > > I know that of course as that is how I loop in my own c or matlab > > programs! One loop over elements and one loop over the faces. > > > I am just saying state in the documentation example > > that plus is on the inside when the normal is pointing outward, > > and minus on the outside. Just add one sentence. > > I usually tend to do that the other way around as > > it is the limit with -epsilon and epsilon positive that is on the inside > (so minus not plus is my usual convention; it does not matter but define it > also > > in the documentation example for clarity). > > > * The one remaining issue in my code seems that FD does not like me to use > D_(+) or D_(-) while D(+) and D(-) are allowed, with the one > > difference being that D is a TrailFunction > > and D_ (previous time step) a Function. > > I guess I don't understand that difference yet. > > > Thanks for the UFL link. That will help. > > > Best wishes, > > > > > > > ------------------------------ > *From:* firedrake-bounces(a)imperial.ac.uk <firedrake-bounces(a)imperial.ac.uk> > on behalf of Colin Cotter <colin.cotter(a)imperial.ac.uk> > *Sent:* Sunday, June 26, 2016 3:11 PM > > *To:* firedrake > *Subject:* Re: [firedrake] upwind example remark; extension to CN in time > > Hi Onno, > This is documented in the UFL manual (the plusses, minusses etc.). > > fenicsproject.org/pub/documents/*ufl*/*ufl*-user-*manual*/*ufl*-user- > *manual*.pdf > > The reason for this notation is that UFL is written for a global domain, > in particular for facets it is the union of all interior facets, and when > surface integrals are done, the facets are visited one by one instead of > the cells. This means that there needs to be an (arbitrary) labelling of + > and - sides of each facet. The notation was originally invented to write > down global inner products to analyse DG methods for elliptic problems, > where it is completely standard. This notation was set out in > Arnold, Douglas N., et al. "Unified analysis of discontinuous Galerkin > methods for elliptic problems." *SIAM journal on numerical analysis* 39.5 > (2002): 1749-1779. > > all the best > --cjc > > On 26 June 2016 at 14:08, Onno Bokhove <O.Bokhove(a)leeds.ac.uk> wrote: > >> Sure I knew that. D+ is on the inside of a cell with n pointing outward >> >> unless I am blinded by being on the other side of the world. >> >> Bit unusual convention in the DG world. Usually minus is on the inside. >> >> But could it as such be defined in the write-up/documentation please? >> >> In any DG article the convention whatever it is gets defined. >> >> >> ------------------------------ >> *From:* firedrake-bounces(a)imperial.ac.uk < >> firedrake-bounces(a)imperial.ac.uk> on behalf of Andrew McRae < >> A.T.T.McRae(a)bath.ac.uk> >> *Sent:* Sunday, June 26, 2016 12:40 PM >> *To:* firedrake(a)imperial.ac.uk >> *Subject:* Re: [firedrake] upwind example remark; extension to CN in time >> >> D+ and D- refer to the values either side of an interior facet between >> two cells. >> >> On the boundary (exterior facets), there is no such ambiguity, so just D >> is used. >> >> In the steady-state example on the website, this D enters the >> left-hand-side integral over the outflow part of the boundary, while the >> prescribed function D0 is used in the integral over the inflow part. >> >> "In general", you can write a term that contains both D and D0 and 'picks >> the right one' depending on the sign of u.n (this is how the interior facet >> integral works). Note that n always points *outwards* on the exterior >> boundary facets. >> >> On 26 June 2016 at 12:24, Onno Bokhove <O.Bokhove(a)leeds.ac.uk> wrote: >> >>> The issue is that D_('+') (or D_('-')) is not defined but I don't >>> understand why not, or how to fix it? >>> >>> >>> ------------------------------ >>> *From:* firedrake-bounces(a)imperial.ac.uk < >>> firedrake-bounces(a)imperial.ac.uk> on behalf of Onno Bokhove < >>> O.Bokhove(a)leeds.ac.uk> >>> *Sent:* Sunday, June 26, 2016 12:21 PM >>> >>> *To:* firedrake >>> *Subject:* Re: [firedrake] upwind example remark; extension to CN in >>> time >>> >>> >>> Yes that was my point indeed. >>> >>> To build any more advanced code upon this one, one >>> >>> would need to use an alternate approach entirely >>> >>> through the flux. >>> >>> That is why it may be useful to extend the >>> >>> example with two alternatives: the simple one shown >>> >>> and one flux-based only. >>> >>> >>> By the way, in the attached I have L2=0 and thet=0; >>> >>> if I uncomment the actual L2 there is a complaint, >>> >>> even though they=0. >>> >>> Obviously I need this term for the general theta-scheme, >>> >>> including thet=0.5 (Crank-Nicolson). >>> >>> I am doing the above simple extensions to get used to FD! >>> Small steps ... >>> >>> The example is ties such that value 1 advects into through the >>> entire domain in the time set. >>> ------------------------------ >>> *From:* firedrake-bounces(a)imperial.ac.uk < >>> firedrake-bounces(a)imperial.ac.uk> on behalf of Colin Cotter < >>> colin.cotter(a)imperial.ac.uk> >>> *Sent:* Sunday, June 26, 2016 11:26 AM >>> *To:* firedrake >>> *Subject:* Re: [firedrake] upwind example remark; extension to CN in >>> time >>> >>> Hi Onno, >>> In answer to your first question. Things are simplified in the >>> advection demo because we know which edges are upwind, which are downwind, >>> and which have no flow through them. If you wanted to generalise this, then >>> we would need a surface integral over the entire external boundary with a >>> switch that either gave the internal value of the tracer or the inflow >>> boundary condition, depending on the direction of u.n, just like we did for >>> the interior facets. >>> Instead, in this example, on the outflow boundary we just write an >>> integral involving the interior value, and on the inflow boundary we write >>> an integral setting the inflow values. >>> >>> all the best >>> -cjc >>> >>> On 25 June 2016 at 16:05, Onno Bokhove <O.Bokhove(a)leeds.ac.uk> wrote: >>> >>>> Hi FDs, >>>> >>>> >>>> (i) The upwind example is odd in that the real boundaries are treated >>>> >>>> in a special way rather than setting the outside value in >>>> >>>> the same numerical flux as anywhere else; >>>> >>>> my question is how this can be fixed? >>>> >>>> So no special terms at the in- and outflow, just specify >>>> >>>> the outside velocity used. This is imo the proper DG way not was is done >>>> >>>> now. >>>> >>>> >>>> (ii) Attached my attempt to do a time dependent advection >>>> >>>> with Crank Nicolson, which works but does not yet >>>> >>>> give the right answer, in part due to (i). >>>> >>>> >>>> [I don't understand why I have to make >>>> >>>> my theta variable that a constant but other wise it did not work.] >>>> >>>> Any suggestions, corrections, also on other silly things done by me? >>>> >>>> Next step is to extend to system with 3 equations; >>>> hence the few extra but unused variables. >>>> >>>> Thanks, from midnight Seoul in old style Korean guesthouse. >>>> >>>> Onno >>>> ------------------------------ >>>> >>> >>> >>> >>> -- >>> http://www.imperial.ac.uk/people/colin.cotter >>> >>> www.cambridge.org/9781107663916 >>> >>> >>> >> > > > -- > http://www.imperial.ac.uk/people/colin.cotter > > www.cambridge.org/9781107663916 > > > -- http://www.imperial.ac.uk/people/colin.cotter www.cambridge.org/9781107663916

2 1

Re: [firedrake] upwind example remark; extension to CN in time
by Colin Cotter 27 Jun '16

27 Jun '16

Hi Onno, This is documented in the UFL manual (the plusses, minusses etc.). fenicsproject.org/pub/documents/*ufl*/*ufl*-user-*manual*/*ufl*-user- *manual*.pdf The reason for this notation is that UFL is written for a global domain, in particular for facets it is the union of all interior facets, and when surface integrals are done, the facets are visited one by one instead of the cells. This means that there needs to be an (arbitrary) labelling of + and - sides of each facet. The notation was originally invented to write down global inner products to analyse DG methods for elliptic problems, where it is completely standard. This notation was set out in Arnold, Douglas N., et al. "Unified analysis of discontinuous Galerkin methods for elliptic problems." *SIAM journal on numerical analysis* 39.5 (2002): 1749-1779. all the best --cjc On 26 June 2016 at 14:08, Onno Bokhove <O.Bokhove(a)leeds.ac.uk> wrote: > Sure I knew that. D+ is on the inside of a cell with n pointing outward > > unless I am blinded by being on the other side of the world. > > Bit unusual convention in the DG world. Usually minus is on the inside. > > But could it as such be defined in the write-up/documentation please? > > In any DG article the convention whatever it is gets defined. > > > ------------------------------ > *From:* firedrake-bounces(a)imperial.ac.uk <firedrake-bounces(a)imperial.ac.uk> > on behalf of Andrew McRae <A.T.T.McRae(a)bath.ac.uk> > *Sent:* Sunday, June 26, 2016 12:40 PM > *To:* firedrake(a)imperial.ac.uk > *Subject:* Re: [firedrake] upwind example remark; extension to CN in time > > D+ and D- refer to the values either side of an interior facet between two > cells. > > On the boundary (exterior facets), there is no such ambiguity, so just D > is used. > > In the steady-state example on the website, this D enters the > left-hand-side integral over the outflow part of the boundary, while the > prescribed function D0 is used in the integral over the inflow part. > > "In general", you can write a term that contains both D and D0 and 'picks > the right one' depending on the sign of u.n (this is how the interior facet > integral works). Note that n always points *outwards* on the exterior > boundary facets. > > On 26 June 2016 at 12:24, Onno Bokhove <O.Bokhove(a)leeds.ac.uk> wrote: > >> The issue is that D_('+') (or D_('-')) is not defined but I don't >> understand why not, or how to fix it? >> >> >> ------------------------------ >> *From:* firedrake-bounces(a)imperial.ac.uk < >> firedrake-bounces(a)imperial.ac.uk> on behalf of Onno Bokhove < >> O.Bokhove(a)leeds.ac.uk> >> *Sent:* Sunday, June 26, 2016 12:21 PM >> >> *To:* firedrake >> *Subject:* Re: [firedrake] upwind example remark; extension to CN in time >> >> >> Yes that was my point indeed. >> >> To build any more advanced code upon this one, one >> >> would need to use an alternate approach entirely >> >> through the flux. >> >> That is why it may be useful to extend the >> >> example with two alternatives: the simple one shown >> >> and one flux-based only. >> >> >> By the way, in the attached I have L2=0 and thet=0; >> >> if I uncomment the actual L2 there is a complaint, >> >> even though they=0. >> >> Obviously I need this term for the general theta-scheme, >> >> including thet=0.5 (Crank-Nicolson). >> >> I am doing the above simple extensions to get used to FD! >> Small steps ... >> >> The example is ties such that value 1 advects into through the >> entire domain in the time set. >> ------------------------------ >> *From:* firedrake-bounces(a)imperial.ac.uk < >> firedrake-bounces(a)imperial.ac.uk> on behalf of Colin Cotter < >> colin.cotter(a)imperial.ac.uk> >> *Sent:* Sunday, June 26, 2016 11:26 AM >> *To:* firedrake >> *Subject:* Re: [firedrake] upwind example remark; extension to CN in time >> >> Hi Onno, >> In answer to your first question. Things are simplified in the >> advection demo because we know which edges are upwind, which are downwind, >> and which have no flow through them. If you wanted to generalise this, then >> we would need a surface integral over the entire external boundary with a >> switch that either gave the internal value of the tracer or the inflow >> boundary condition, depending on the direction of u.n, just like we did for >> the interior facets. >> Instead, in this example, on the outflow boundary we just write an >> integral involving the interior value, and on the inflow boundary we write >> an integral setting the inflow values. >> >> all the best >> -cjc >> >> On 25 June 2016 at 16:05, Onno Bokhove <O.Bokhove(a)leeds.ac.uk> wrote: >> >>> Hi FDs, >>> >>> >>> (i) The upwind example is odd in that the real boundaries are treated >>> >>> in a special way rather than setting the outside value in >>> >>> the same numerical flux as anywhere else; >>> >>> my question is how this can be fixed? >>> >>> So no special terms at the in- and outflow, just specify >>> >>> the outside velocity used. This is imo the proper DG way not was is done >>> >>> now. >>> >>> >>> (ii) Attached my attempt to do a time dependent advection >>> >>> with Crank Nicolson, which works but does not yet >>> >>> give the right answer, in part due to (i). >>> >>> >>> [I don't understand why I have to make >>> >>> my theta variable that a constant but other wise it did not work.] >>> >>> Any suggestions, corrections, also on other silly things done by me? >>> >>> Next step is to extend to system with 3 equations; >>> hence the few extra but unused variables. >>> >>> Thanks, from midnight Seoul in old style Korean guesthouse. >>> >>> Onno >>> ------------------------------ >>> >> >> >> >> -- >> http://www.imperial.ac.uk/people/colin.cotter >> >> www.cambridge.org/9781107663916 >> >> >> > -- http://www.imperial.ac.uk/people/colin.cotter www.cambridge.org/9781107663916

2 1