- firedrake - imperial.ac.uk

Re: [firedrake] Mixed system
by Colin Cotter 12 Oct '16

12 Oct '16

OK, I will try to discuss this with Onno on Thursday once he has finished being Examiner In Chief. cheers --cjc On 10 October 2016 at 18:29, Mitchell, Lawrence < lawrence.mitchell(a)imperial.ac.uk> wrote: > > > On 10 Oct 2016, at 18:19, Colin Cotter <colin.cotter(a)imperial.ac.uk> > wrote: > > > > OK, this is the problem then - you need to find an appropriate > preconditioner. > > > > You need to start by finding an appropriate preconditioner for the > system without the R space part. > > > > Lawrence: once that is done, is there anything special to do to > incorporate it into the system with the R space part [asking because you > seem to have thought about this previously], or do we just recommend block > preconditioning in this case? > > I think you could use block first. It's easier. Otherwise I think the > firedrake webpages sketch how to do sherman-morrison. > > Lawrence > _______________________________________________ > 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 2

Re: [firedrake] Mixed system
by Colin Cotter 07 Oct '16

07 Oct '16

We probably need to think about the preconditioner. First please could you confirm that the timing is dominated by the linear solver? On Wednesday, 5 October 2016, Onno Bokhove <O.Bokhove(a)leeds.ac.uk> wrote: > Great news! Curious about the next step ...!?! > > > ------------------------------ > *From:* firedrake-bounces(a)imperial.ac.uk > <javascript:_e(%7B%7D,'cvml','firedrake-bounces(a)imperial.ac.uk');> < > firedrake-bounces(a)imperial.ac.uk > <javascript:_e(%7B%7D,'cvml','firedrake-bounces(a)imperial.ac.uk');>> on > behalf of Anna Kalogirou <a.kalogirou(a)leeds.ac.uk > <javascript:_e(%7B%7D,'cvml','a.kalogirou(a)leeds.ac.uk');>> > *Sent:* Wednesday, October 5, 2016 6:43 PM > *To:* firedrake(a)imperial.ac.uk > <javascript:_e(%7B%7D,'cvml','firedrake(a)imperial.ac.uk');> > *Subject:* Re: [firedrake] Mixed system > > Dear all, > > So now that the real function spaces are up and running (thanks for this!) > I can confirm that my code runs without errors and gives the same result as > the one I had before, which was using linear algebra techniques to take > care of the constraint enforced with a Lagrange multiplier. > > However, I noticed that the old code > <https://bitbucket.org/annakalog/buoy2d/src> (without mixed function > spaces) takes 27.5 minutes for a given final time T, while the new > <https://bitbucket.org/annakalog/buoy2d/src/a1760633fb44f99b3b5ba4d09682b7c7…> > code takes almost twice as much to run up to the same T, namely 47.5 > minutes. > > Any ideas why this is happening? > > Best wishes, > > Anna. > > > On 16/09/16 17:22, David Ham wrote: > > Hi Anna, > > Are you certain that the test you are running is the one you pointed to on > the web. Line 89 from your error message does not appear in the version on > bitbucket, and the last commit was 4 days ago. > > Regards, > > David > > On Fri, 16 Sep 2016 at 17:03 Anna Kalogirou <A.Kalogirou(a)leeds.ac.uk > <javascript:_e(%7B%7D,'cvml','A.Kalogirou(a)leeds.ac.uk');>> wrote: > >> When I run >> >> firedrake-zenodo -m "Versions I am using now” >> >> I get >> >> Failed to retrieve git hash for optional component 'SLOPE', continuing >> without it >> >> The file firedrake.json contains the following (up to ufl it’s the same >> as yours): >> >> {"fiat": "1bc440ca34d4b127897cf5a36d63bf8a41a3d55c", "COFFEE": " >> 8e03328d659fd8d8871bf5b7bc3f08b769f0ab9d", "ufl": " >> 2fc6d00f7e4b5e0639c6d01f27f9b2bbc0d4696e", "firedrake": " >> 71cf933f05508427287dc6464e15d2d7adf1b179", "tsfc": " >> f867938c58f475f1cd332421e2bf14cffa2f5aef", "petsc4py": " >> 0b590afa9d3c42e0b4f863061b6e082d6f0227f6", "petsc": " >> 21143a8c817f512a23f8e18da648891eb40f5217", "PyOP2": " >> d9ea13e1e17fcb1e05adcf1db31bd8190974520d", "message": "Versions I am >> using now"} >> >> >> On 16 Sep 2016, at 16:49, David Ham <David.Ham(a)imperial.ac.uk >> <javascript:_e(%7B%7D,'cvml','David.Ham(a)imperial.ac.uk');>> wrote: >> >> OK. Let's do some sanity checking. >> >> 1. Please confirm that the current firedrake checkout is in developer >> mode and has always been in developer mode (I can see that it is in >> developer mode from the error message but I want to make sure there cannot >> be stale files around from a non-developer mode past). >> >> 2. Run firedrake-clean and check that the error still occurs. >> >> 2. Please run the following command: >> >> firedrake-zenodo -m "Versions I am using now" >> >> This will create a file firedrake.json . Please check that its contents >> exactly match what is below: >> >> {"fiat": "1bc440ca34d4b127897cf5a36d63bf8a41a3d55c", "COFFEE": " >> 8e03328d659fd8d8871bf5b7bc3f08b769f0ab9d", "ufl": " >> 2fc6d00f7e4b5e0639c6d01f27f9b2bbc0d4696e", "firedrake": " >> 71cf933f05508427287dc6464e15d2d7adf1b179", "tsfc": " >> f867938c58f475f1cd332421e2bf14cffa2f5aef", "petsc4py": " >> 0b590afa9d3c42e0b4f863061b6e082d6f0227f6", "petsc": " >> 21143a8c817f512a23f8e18da648891eb40f5217", "PyOP2": " >> d9ea13e1e17fcb1e05adcf1db31bd8190974520d", "message": "Versions I am >> using now"} >> >> >> On Fri, 16 Sep 2016 at 16:28 Anna Kalogirou <A.Kalogirou(a)leeds.ac.uk >> <javascript:_e(%7B%7D,'cvml','A.Kalogirou(a)leeds.ac.uk');>> wrote: >> >>> That is the same version I have. >>> >>> >>> On 16 Sep 2016, at 16:23, David Ham <David.Ham(a)imperial.ac.uk >>> <javascript:_e(%7B%7D,'cvml','David.Ham(a)imperial.ac.uk');>> wrote: >>> >>> Hi Anna, >>> >>> Please confirm the git revision of your PyOP2. You can retrieve this by >>> typing the following command in the PyOP2 directory: >>> >>> git rev-parse HEAD >>> >>> You should have: >>> d9ea13e1e17fcb1e05adcf1db31bd8190974520d >>> >>> Regards, >>> >>> David >>> >>> On Fri, 16 Sep 2016 at 15:59 Anna Kalogirou <A.Kalogirou(a)leeds.ac.uk >>> <javascript:_e(%7B%7D,'cvml','A.Kalogirou(a)leeds.ac.uk');>> wrote: >>> >>>> Hi David, >>>> >>>> Thank you for the reply. Indeed, the problem was with the virtualenv; I >>>> had to deactivate and activate it again. Unfortunately I still get errors, >>>> which are probably due to the firedrake I have installed (since you said >>>> that it works on your computer): >>>> >>>> File "buoy-swe.py", line 89, in <module> >>>> F_solver = solver_F(phi0_5, eta1, mu0_5, I, w, phi0, eta0, Z0, W0, >>>> step_b, etaR, phi_t, eta_t, mu_t, I_t, v1, v2, v3, v4, dt, Hb, H0, L, >>>> dR_dt, g, rho, Mass, solvers_print); >>>> File "/Users/matak/Documents/Simulations/Firedrake/Ship/Modules/Mixed >>>> system/solvers.py", line 41, in solver_F >>>> F_solver = LinearVariationalSolver(F_problem, >>>> solver_parameters=solvers_print) >>>> File "/Users/matak/firedrake/src/firedrake/firedrake/variational_solver.py", >>>> line 264, in __init__ >>>> super(LinearVariationalSolver, self).__init__(*args, **kwargs) >>>> File "/Users/matak/firedrake/src/firedrake/firedrake/variational_solver.py", >>>> line 127, in __init__ >>>> ctx = solving_utils._SNESContext(problem) >>>> File "/Users/matak/firedrake/src/firedrake/firedrake/solving_utils.py", >>>> line 107, in __init__ >>>> for problem in problems) >>>> File "/Users/matak/firedrake/src/firedrake/firedrake/solving_utils.py", >>>> line 107, in <genexpr> >>>> for problem in problems) >>>> File "/Users/matak/firedrake/src/firedrake/firedrake/assemble.py", >>>> line 67, in assemble >>>> inverse=inverse, nest=nest) >>>> File "<decorator-gen-294>", line 2, in _assemble >>>> File "/Users/matak/firedrake/src/firedrake/firedrake/utils.py", line >>>> 62, in wrapper >>>> return f(*args, **kwargs) >>>> File "/Users/matak/firedrake/src/firedrake/firedrake/assemble.py", >>>> line 179, in _assemble >>>> nest=nest) >>>> File "/Users/matak/firedrake/src/PyOP2/pyop2/backends.py", line 118, >>>> in __call__ >>>> return t(*args, **kwargs) >>>> File "/Users/matak/firedrake/src/PyOP2/pyop2/caching.py", line 122, >>>> in __new__ >>>> args, kwargs = cls._process_args(*args, **kwargs) >>>> File "<decorator-gen-258>", line 2, in _process_args >>>> File "/Users/matak/firedrake/src/PyOP2/pyop2/utils.py", line 130, in >>>> wrapper >>>> return f(*args, **kwargs) >>>> File "/Users/matak/firedrake/src/PyOP2/pyop2/base.py", line 3591, in >>>> _process_args >>>> if not (pair[0].toset == dsets[0].set and >>>> File "/Users/matak/firedrake/src/PyOP2/pyop2/utils.py", line 64, in >>>> __get__ >>>> obj.__dict__[self.__name__] = result = self.fget(obj) >>>> File "/Users/matak/firedrake/src/PyOP2/pyop2/base.py", line 3372, in >>>> toset >>>> m.toset for m in self._maps)) >>>> File "/Users/matak/firedrake/src/PyOP2/pyop2/caching.py", line 156, >>>> in __new__ >>>> obj = make_obj() >>>> File "/Users/matak/firedrake/src/PyOP2/pyop2/caching.py", line 137, >>>> in make_obj >>>> obj.__init__(*args, **kwargs) >>>> File "/Users/matak/firedrake/src/PyOP2/pyop2/base.py", line 977, in >>>> __init__ >>>> "All components of a MixedSet must have the same number of layers." >>>> AssertionError: All components of a MixedSet must have the same number >>>> of layers. >>>> >>>> >>>> >>>> On 16 Sep 2016, at 15:47, David Ham <David.Ham(a)imperial.ac.uk >>>> <javascript:_e(%7B%7D,'cvml','David.Ham(a)imperial.ac.uk');>> wrote: >>>> >>>> Hi Anna, >>>> >>>> Your code works for me (in the sense that it appears to run without >>>> error. I have no idea whether the answer is right ;) >>>> >>>> Please also note that Miklos is correct: your error is definitely from >>>> a global firedrake install. You can tell because the paths start >>>> /usr/local. Your previous error messages are from a local install. You can >>>> tell that because the paths start /Users/matak . Are you sure your >>>> virtualenv is active? >>>> >>>> Regards, >>>> >>>> David >>>> >>>> >>>> >>>> On Fri, 16 Sep 2016 at 15:16 Anna Kalogirou <A.Kalogirou(a)leeds.ac.uk >>>> <javascript:_e(%7B%7D,'cvml','A.Kalogirou(a)leeds.ac.uk');>> wrote: >>>> >>>>> During the past few days I have been working on my laptop (not the >>>>> desktop at Leeds) and I don’t believe there are multiple installations of >>>>> firedrake on this machine. >>>>> >>>>> Can someone please try to run my code >>>>> <https://bitbucket.org/annakalog/buoy2d/src/673fef3466b06c98ac11f5a120ff4211…> so >>>>> that we know for sure if the problem is in my code or related with the >>>>> installation? Thanks. >>>>> >>>>> >>>>> On 16 Sep 2016, at 14:54, Colin Cotter <colin.cotter(a)imperial.ac.uk >>>>> <javascript:_e(%7B%7D,'cvml','colin.cotter(a)imperial.ac.uk');>> wrote: >>>>> >>>>> Ah yes, that tallies with Daniel Ruprecht's reports that Firedrake is >>>>> part of the standard installation at Leeds! >>>>> >>>>> On 16 September 2016 at 14:53, Homolya, Miklós < >>>>> m.homolya14(a)imperial.ac.uk> wrote: >>>>> >>>>> It seems that the Firedrake you're importing is a system wide >>>>> installation. No wonder that updating your own installation doesn't fix it. >>>>> >>>>> Get Outlook for Android <https://aka.ms/ghei36> >>>>> >>>>> >>>>> >>>>> >>>>> On Thu, Sep 15, 2016 at 8:33 PM +0100, "Anna Kalogirou" < >>>>> A.Kalogirou(a)leeds.ac.uk> wrote: >>>>> >>>>> Ok thanks, I get the error below: >>>>> >>>>> Traceback (most recent call last): >>>>> >>>>> File "buoy-swe.py", line 6, in <module> >>>>> from firedrake import * >>>>> File "/usr/local/lib/python2.7/site-packages/firedrake/__init__.py", >>>>> line 30, in <module> >>>>> from assemble import * >>>>> File "/usr/local/lib/python2.7/site-packages/firedrake/assemble.py", >>>>> line 8, in <module> >>>>> import assembly_cache >>>>> File "/usr/local/lib/python2.7/site-packages/firedrake/assembly_cache.py", >>>>> line 40, in <module> >>>>> import function >>>>> File "/usr/local/lib/python2.7/site-packages/firedrake/function.py", >>>>> line 9, in <module> >>>>> import assemble_expressions >>>>> File "/usr/local/lib/python2.7/site-packages/firedrake/assemble_expressions.py", >>>>> line 15, in <module> >>>>> import functionspace >>>>> File "/usr/local/lib/python2.7/site-packages/firedrake/functionspace.py", >>>>> line 10, in <module> >>>>> import dmplex >>>>> ImportError: dlopen(/usr/local/lib/python2.7/site-packages/firedrake/dmplex.so, >>>>> 2): Library not loaded: /usr/local/lib/python2.7/site- >>>>> packages/petsc/lib/libpetsc.3.05.2.dylib >>>>> Referenced from: /usr/local/lib/python2.7/site- >>>>> packages/firedrake/dmplex.so >>>>> Reason: image not found >>>>> >>>>> On 15 Sep 2016, at 18:21, Andrew McRae <A.T.T.McRae(a)bath.ac.uk> wrote: >>>>> >>>>> To give credit where it's due, >>>>> >>>>> s/Andrew's patch/Miklós's patch/ >>>>> >>>>> On 15 September 2016 at 18:16, David Ham <David.Ham(a)imperial.ac.uk> >>>>> wrote: >>>>> >>>>> I suspect Andrew is correct. However the deeper issue is that you >>>>> edited the installed file, which is never a safe thing to do. One must >>>>> always edit the file in the source directory and then (if necessary) >>>>> re-install. >>>>> >>>>> Anyway, I have committed Andrew's patch to the repository, so >>>>> hopefully the problem will go away if you do a git update. >>>>> >>>>> Regards, >>>>> >>>>> David >>>>> >>>>> On Wed, 14 Sep 2016 at 18:39 Andrew McRae <A.T.T.McRae(a)bath.ac.uk> >>>>> wrote: >>>>> >>>>> Probably related to installing Firedrake in development mode or not. >>>>> >>>>> I can't help any further here; I don't know if someone else can. >>>>> >>>>> On 14 September 2016 at 16:13, Anna Kalogirou <A.Kalogirou(a)leeds.ac.uk >>>>> > wrote: >>>>> >>>>> Ok, so I replace line 3366 in file “/Users/matak/firedrake/lib/ >>>>> python2.7/site-packages/pyop2/base.py” with the suggested patch, but >>>>> I now get the error below. >>>>> >>>>> Note that the error I get comes from file firedrake/lib/python2.7/ >>>>> site-packages/pyop2/base.py >>>>> while in the reported issue the error comes from firedrake/src/PyOP2/ >>>>> pyop2/base.py >>>>> so the two might not be related after all. >>>>> >>>>> Thanks, >>>>> Anna. >>>>> >>>>> >>>>> >>>>> Traceback (most recent call last): >>>>> File "buoy-swe.py", line 89, in <module> >>>>> F_solver = solver_F(phi0_5, eta1, mu0_5, I, w, phi0, eta0, Z0, W0, >>>>> step_b, etaR, phi_t, eta_t, mu_t, I_t, v1, v2, v3, v4, dt, Hb, H0, L, >>>>> dR_dt, g, rho, Mass, solvers_print); >>>>> File "/Users/matak/Documents/Simulations/Firedrake/Ship/Modules/Mixed >>>>> system/solvers.py", line 41, in solver_F >>>>> F_solver = LinearVariationalSolver(F_problem, >>>>> solver_parameters=solvers_print) >>>>> File "/Users/matak/firedrake/lib/python2.7/site-packages/ >>>>> firedrake/variational_solver.py", line 264, in __init__ >>>>> super(LinearVariationalSolver, self).__init__(*args, **kwargs) >>>>> File "/Users/matak/firedrake/lib/python2.7/site-packages/ >>>>> firedrake/variational_solver.py", line 127, in __init__ >>>>> ctx = solving_utils._SNESContext(problem) >>>>> File "/Users/matak/firedrake/lib/python2.7/site-packages/firedrake/solving_utils.py", >>>>> line 107, in __init__ >>>>> for problem in problems) >>>>> File "/Users/matak/firedrake/lib/python2.7/site-packages/firedrake/solving_utils.py", >>>>> line 107, in <genexpr> >>>>> for problem in problems) >>>>> File "/Users/matak/firedrake/lib/python2.7/site-packages/firedrake/assemble.py", >>>>> line 67, in assemble >>>>> inverse=inverse, nest=nest) >>>>> File "<decorator-gen-294>", line 2, in _assemble >>>>> File "/Users/matak/firedrake/lib/python2.7/site-packages/firedrake/utils.py", >>>>> line 62, in wrapper >>>>> return f(*args, **kwargs) >>>>> File "/Users/matak/firedrake/lib/python2.7/site-packages/firedrake/assemble.py", >>>>> line 179, in _assemble >>>>> nest=nest) >>>>> File "/Users/matak/firedrake/lib/python2.7/site-packages/pyop2/backends.py", >>>>> line 118, in __call__ >>>>> return t(*args, **kwargs) >>>>> File "/Users/matak/firedrake/lib/python2.7/site-packages/pyop2/caching.py", >>>>> line 122, in __new__ >>>>> args, kwargs = cls._process_args(*args, **kwargs) >>>>> File "<decorator-gen-258>", line 2, in _process_args >>>>> File "/Users/matak/firedrake/lib/python2.7/site-packages/pyop2/utils.py", >>>>> line 130, in wrapper >>>>> return f(*args, **kwargs) >>>>> File "/Users/matak/firedrake/lib/python2.7/site-packages/pyop2/base.py", >>>>> line 3592, in _process_args >>>>> if not (pair[0].toset == dsets[0].set and >>>>> File "/Users/matak/firedrake/lib/python2.7/site-packages/pyop2/utils.py", >>>>> line 64, in __get__ >>>>> obj.__dict__[self.__name__] = result = self.fget(obj) >>>>> File "/Users/matak/firedrake/lib/python2.7/site-packages/pyop2/base.py", >>>>> line 3373, in toset >>>>> m.toset for m in self._maps)) >>>>> File "/Users/matak/firedrake/lib/python2.7/site-packages/pyop2/caching.py", >>>>> line 156, in __new__ >>>>> obj = make_obj() >>>>> File "/Users/matak/firedrake/lib/python2.7/site-packages/pyop2/caching.py", >>>>> line 137, in make_obj >>>>> obj.__init__(*args, **kwargs) >>>>> File "/Users/matak/firedrake/lib/python2.7/site-packages/pyop2/base.py", >>>>> line 977, in __init__ >>>>> "All components of a MixedSet must have the same number of layers." >>>>> AssertionError: All components of a MixedSet must have the same number >>>>> of layers. >>>>> >>>>> >>>>> ********************************************************** >>>>> >>>>> >>>>> Dr Anna Kalogirou >>>>> Research Fellow >>>>> >>>>> Department of Applied Mathematics >>>>> University of Leeds >>>>> >>>>> Web: http://www1.maths.leeds.ac.uk/~matak/ >>>>> >>>>> Youtube: https://www.youtube.com/channel/UCUZAYHtVoiMqepQflisp66g >>>>> >>>>> ********************************************************** >>>>> Leverhulme Trust Early Career Fellow 2017-2019 >>>>> >>>>> ********************************************************** >>>>> >>>>> >>>>> On 14 Sep 2016, at 16:00, Andrew McRae <A.T.T.McRae(a)bath.ac.uk> wrote: >>>>> >>>>> Not as a branch, but you can change PyOP2 locally using the code at >>>>> the bottom of that post, "Suggested patch" >>>>> >>>>> On 14 September 2016 at 15:54, Anna Kalogirou <A.Kalogirou(a)leeds.ac.uk >>>>> > wrote: >>>>> >>>>> Hi Andrew, >>>>> >>>>> Yes indeed, the two might be related. Was there a fix for the issue >>>>> reported on github? >>>>> >>>>> Best, Anna. >>>>> >>>>> >>>>> ********************************************************** >>>>> >>>>> >>>>> Dr Anna Kalogirou >>>>> Research Fellow >>>>> >>>>> Department of Applied Mathematics >>>>> University of Leeds >>>>> >>>>> Web: http://www1.maths.leeds.ac.uk/~matak/ >>>>> >>>>> Youtube: https://www.youtube.com/channel/UCUZAYHtVoiMqepQflisp66g >>>>> >>>>> ********************************************************** >>>>> Leverhulme Trust Early Career Fellow 2017-2019 >>>>> >>>>> ********************************************************** >>>>> >>>>> >>>>> On 13 Sep 2016, at 16:34, Andrew McRae <A.T.T.McRae(a)bath.ac.uk> wrote: >>>>> >>>>> Probably related to https://github.com/firedrakeproject/firedrake/ >>>>> issues/862 (third post, not first). You might be able to adapt the >>>>> fix; I haven't looked in detail. >>>>> >>>>> On 12 September 2016 at 16:34, Anna Kalogirou <A.Kalogirou(a)leeds.ac.uk >>>>> > wrote: >>>>> >>>>> Dear all, >>>>> >>>>> Regarding the problem described in my email sent of Friday, I fixed a >>>>> minor typo and now get the error below. Any ideas what might be causing it? >>>>> Thanks! >>>>> >>>>> Best, Anna. >>>>> >>>>> Traceback (most recent call last): >>>>> File "buoy-swe.py", line 89, in <module> >>>>> F_solver = solver_F(phi0_5, eta1, mu0_5, I, w, phi0, eta0, Z0, W0, >>>>> step_b, etaR, phi_t, eta_t, mu_t, I_t, v1, v2, v3, v4, dt, Hb, H0, L, >>>>> dR_dt, g, rho, Mass, solvers_print); >>>>> File "/Users/matak/Documents/Simulations/Firedrake/Ship/Modules/Mixed >>>>> system/solvers.py", line 41, in solver_F >>>>> F_solver = LinearVariationalSolver(F_problem, >>>>> solver_parameters=solvers_print) >>>>> File "/Users/matak/firedrake/lib/python2.7/site-packages/ >>>>> firedrake/variational_solver.py", line 264, in __init__ >>>>> super(LinearVariationalSolver, self).__init__(*args, **kwargs) >>>>> File "/Users/matak/firedrake/lib/python2.7/site-packages/ >>>>> firedrake/variational_solver.py", line 127, in __init__ >>>>> ctx = solving_utils._SNESContext(problem) >>>>> File "/Users/matak/firedrake/lib/python2.7/site-packages/firedrake/solving_utils.py", >>>>> line 107, in __init__ >>>>> for problem in problems) >>>>> File "/Users/matak/firedrake/lib/python2.7/site-packages/firedrake/solving_utils.py", >>>>> line 107, in <genexpr> >>>>> for problem in problems) >>>>> File "/Users/matak/firedrake/lib/python2.7/site-packages/firedrake/assemble.py", >>>>> line 67, in assemble >>>>> inverse=inverse, nest=nest) >>>>> File "<decorator-gen-294>", line 2, in _assemble >>>>> File "/Users/matak/firedrake/lib/python2.7/site-packages/firedrake/utils.py", >>>>> line 62, in wrapper >>>>> return f(*args, **kwargs) >>>>> File "/Users/matak/firedrake/lib/python2.7/site-packages/firedrake/assemble.py", >>>>> line 163, in _assemble >>>>> map_pairs.append((op2.DecoratedMap(test.cell_node_map(), >>>>> cell_domains), >>>>> File "/Users/matak/firedrake/lib/python2.7/site-packages/ >>>>> firedrake/functionspaceimpl.py", line 646, in cell_node_map >>>>> for i, s in enumerate(self._spaces)) >>>>> File "/Users/matak/firedrake/lib/python2.7/site-packages/pyop2/backends.py", >>>>> line 118, in __call__ >>>>> return t(*args, **kwargs) >>>>> File "/Users/matak/firedrake/lib/python2.7/site-packages/pyop2/caching.py", >>>>> line 156, in __new__ >>>>> obj = make_obj() >>>>> File "/Users/matak/firedrake/lib/python2.7/site-packages/pyop2/caching.py", >>>>> line 137, in make_obj >>>>> obj.__init__(*args, **kwargs) >>>>> File "/Users/matak/firedrake/lib/python2.7/site-packages/pyop2/base.py", >>>>> line 3342, in __init__ >>>>> if not all(m is None or m.iterset == self.iterset for m in >>>>> self._maps): >>>>> File "/Users/matak/firedrake/lib/python2.7/site-packages/pyop2/base.py", >>>>> line 3342, in <genexpr> >>>>> if not all(m is None or m.iterset == self.iterset for m in >>>>> self._maps): >>>>> File "/Users/matak/firedrake/lib/python2.7/site-packages/pyop2/utils.py", >>>>> line 64, in __get__ >>>>> obj.__dict__[self.__name__] = result = self.fget(obj) >>>>> File "/Users/matak/firedrake/lib/python2.7/site-packages/pyop2/base.py", >>>>> line 3366, in iterset >>>>> return reduce(lambda a, b: a if a is None else a.iterset or b if b >>>>> is None else b.iterset, self._maps) >>>>> File "/Users/matak/firedrake/lib/python2.7/site-packages/pyop2/base.py", >>>>> line 3366, in <lambda> >>>>> return reduce(lambda a, b: a if a is None else a.iterset or b if b >>>>> is None else b.iterset, self._maps) >>>>> File "/Users/matak/firedrake/lib/python2.7/site-packages/pyop2/base.py", >>>>> line 799, in __getattr__ >>>>> return getattr(self._parent, name) >>>>> AttributeError: 'Set' object has no attribute ‘iterset' >>>>> >>>>> >>>>> ********************************************************** >>>>> >>>>> >>>>> Dr Anna Kalogirou >>>>> Research Fellow >>>>> >>>>> Department of Applied Mathematics >>>>> University of Leeds >>>>> >>>>> Web: http://www1.maths.leeds.ac.uk/~matak/ >>>>> >>>>> Youtube: https://www.youtube.com/channel/UCUZAYHtVoiMqepQflisp66g >>>>> >>>>> ********************************************************** >>>>> Leverhulme Trust Early Career Fellow 2017-2019 >>>>> >>>>> ********************************************************** >>>>> >>>>> >>>>> On 9 Sep 2016, at 17:37, Anna Kalogirou <a.kalogirou(a)leeds.ac.uk> >>>>> wrote: >>>>> >>>>> Hi David, >>>>> >>>>> I have tried to implement my problem (attached) using the real >>>>> function spaces but I get the error below. The code can be found here >>>>> <https://bitbucket.org/annakalog/buoy2d/src/d6a805b6a8d0/Mixed%20system/?at=…> >>>>> . >>>>> >>>>> I define the real space by V_R = FunctionSpace(mesh, "R", 0). Is that >>>>> correct? >>>>> Also, is it a problem that the linear part of the variational form >>>>> does not depend on one of the test functions from the mixed function space >>>>> (only depends on v1, v2, v3 and not v4)? >>>>> >>>>> Best, >>>>> >>>>> Anna. >>>>> >>>>> >>>>> >>>>> Traceback (most recent call last): >>>>> File "buoy-swe.py", line 89, in <module> >>>>> F_solver = solver_F(phi0_5, eta1, mu0_5, I, w, phi0, eta0, Z0, W0, >>>>> step_b, etaR, phi_t, eta_t, mu_t, I_t, v1, v2, v3, v4, dt, Hb, H0, L, >>>>> dR_dt, g, rho, Mass, solvers_print); >>>>> File "/Users/matak/Documents/Simulations/Firedrake/Ship/Modules/Mixed >>>>> system/solvers.py", line 41, in solver_F >>>>> F_solver = LinearVariationalSolver(F_problem, >>>>> solver_parameters=solvers_print) >>>>> File "/Users/matak/firedrake/lib/python2.7/site-packages/ >>>>> firedrake/variational_solver.py", line 264, in __init__ >>>>> super(LinearVariationalSolver, self).__init__(*args, **kwargs) >>>>> File "/Users/matak/firedrake/lib/python2.7/site-packages/ >>>>> firedrake/variational_solver.py", line 127, in __init__ >>>>> ctx = solving_utils._SNESContext(problem) >>>>> File "/Users/matak/firedrake/lib/python2.7/site-packages/firedrake/solving_utils.py", >>>>> line 107, in __init__ >>>>> for problem in problems) >>>>> File "/Users/matak/firedrake/lib/python2.7/site-packages/firedrake/solving_utils.py", >>>>> line 107, in <genexpr> >>>>> for problem in problems) >>>>> File "/Users/matak/firedrake/lib/python2.7/site-packages/firedrake/assemble.py", >>>>> line 67, in assemble >>>>> inverse=inverse, nest=nest) >>>>> File "<decorator-gen-294>", line 2, in _assemble >>>>> File "/Users/matak/firedrake/lib/python2.7/site-packages/firedrake/utils.py", >>>>> line 62, in wrapper >>>>> return f(*args, **kwargs) >>>>> File "/Users/matak/firedrake/lib/python2.7/site-packages/firedrake/assemble.py", >>>>> line 101, in _assemble >>>>> inverse=inverse) >>>>> File "/Users/matak/firedrake/lib/python2.7/site-packages/firedrake/tsfc_interface.py", >>>>> line 307, in compile_form >>>>> number_map).kernels >>>>> File "/Users/matak/firedrake/lib/python2.7/site-packages/pyop2/caching.py", >>>>> line 198, in __new__ >>>>> obj = make_obj() >>>>> File "/Users/matak/firedrake/lib/python2.7/site-packages/pyop2/caching.py", >>>>> line 188, in make_obj >>>>> obj.__init__(*args, **kwargs) >>>>> File "/Users/matak/firedrake/lib/python2.7/site-packages/firedrake/tsfc_interface.py", >>>>> line 212, in __init__ >>>>> tree = tsfc_compile_form(form, prefix=name, parameters=parameters) >>>>> File "/Users/matak/firedrake/lib/python2.7/site-packages/tsfc/driver.py", >>>>> line 47, in compile_form >>>>> do_estimate_degrees=True) >>>>> File "/Users/matak/firedrake/lib/python2.7/site-packages/ufl/ >>>>> algorithms/compute_form_data.py", line 386, in compute_form_data >>>>> check_form_arity(preprocessed_form, self.original_form.arguments()) >>>>> # Currently testing how fast this is >>>>> File "/Users/matak/firedrake/lib/python2.7/site-packages/ufl/algorithms/check_arities.py", >>>>> line 152, in check_form_arity >>>>> check_integrand_arity(itg.integrand(), arguments) >>>>> File "/Users/matak/firedrake/lib/python2.7/site-packages/ufl/algorithms/check_arities.py", >>>>> line 145, in check_integrand_arity >>>>> args = map_expr_dag(rules, expr, compress=False) >>>>> File "/Users/matak/firedrake/lib/python2.7/site-packages/ufl/corealg/map_dag.py", >>>>> line 37, in map_expr_dag >>>>> result, = map_expr_dags(function, [expression], compress=compress) >>>>> File "/Users/matak/firedrake/lib/python2.7/site-packages/ufl/corealg/map_dag.py", >>>>> line 86, in map_expr_dags >>>>> r = handlers[v._ufl_typecode_](v, *[vcache[u] for u in >>>>> v.ufl_operands]) >>>>> File "/Users/matak/firedrake/lib/python2.7/site-packages/ufl/algorithms/check_arities.py", >>>>> line 42, in sum >>>>> raise ArityMismatch("Adding expressions with non-matching form >>>>> arguments {0} vs {1}.".format(a, b)) >>>>> ufl.algorithms.check_arities.ArityMismatch: Adding expressions with >>>>> non-matching form arguments () vs (Argument(WithGeometry( >>>>> IndexedProxyFunctionSpace(<firedrake.mesh.ExtrudedMeshTopology object >>>>> at 0x113e6b890>, TensorProductElement(FiniteElement('Lagrange', >>>>> interval, 1), FiniteElement('Lagrange', interval, 1), >>>>> cell=TensorProductCell(interval, interval)), name=None, index=0, >>>>> component=None), Mesh(VectorElement(TensorProductElement(FiniteElement('Lagrange', >>>>> interval, 1), FiniteElement('Lagrange', interval, 1), >>>>> cell=TensorProductCell(interval, interval)), dim=2), 4)), 1, None),). >>>>> <coupled_system.pdf> >>>>> >>>>> On 18 Aug 2016, at 16:03, David Ham <david.ham(a)imperial.ac.uk> wrote: >>>>> >>>>> Hi Anna, >>>>> >>>>> Real function spaces now pass some simple tests, at least on my >>>>> machine. You need the real-function-space branch of firedrake, and the >>>>> globals_matrices branch of PyOP2. It's not going to get merged before I >>>>> disappear for a week - not least because I need to update documentation. >>>>> However you should be able to test it out on the branches. >>>>> >>>>> Regards, >>>>> >>>>> David >>>>> >>>>> On Wed, 10 Aug 2016 at 16:58 David Ham <David.Ham(a)imperial.ac.uk> >>>>> wrote: >>>>> >>>>> Unfortunately, It isn't as simple as saying "no mpi". MPI pervades >>>>> Firedrake. Unless real function spaces play nicely with communicators, >>>>> Firedrake will simply crash. I'm sorry that this feature is not available >>>>> right now, but it actually does require significant work. >>>>> >>>>> On Wed, 10 Aug 2016 at 16:51 Anna Kalogirou <A.Kalogirou(a)leeds.ac.uk> >>>>> wrote: >>>>> >>>>> Yes indeed, but I believe David would want to make it work for all >>>>> cases. The merge is necessary in order to proceed with my simulations, so >>>>> an update as soon as possible would be greatly appreciated. >>>>> >>>>> Best, Anna. >>>>> >>>>> >>>>> On 10 Aug 2016, at 16:29, Onno Bokhove <O.Bokhove(a)leeds.ac.uk> wrote: >>>>> >>>>> We would already be helped for the case without MPI, I think that is >>>>> correct, Anna? >>>>> >>>>> >>>>> *From:* firedrake-bounces(a)imperial.ac.uk <firedrake-bounces@imperial. >>>>> ac.uk> on behalf of David Ham <David.Ham(a)imperial.ac.uk> >>>>> >>>>> *Sent:* Wednesday, August 10, 2016 10:01 AM >>>>> >>>>> >>>>> *To:* firedrake >>>>> *Subject:* Re: [firedrake] Mixed system >>>>> >>>>> The short version is it depends on when I find some time to work on >>>>> it. It is not simply the case that we have new features lying around >>>>> waiting to be merged. The current proximate issue is to make the real >>>>> spaces work properly with MPI communicators. I can't guarantee that will be >>>>> the last issue. >>>>> >>>>> >>>>> On Tue, 9 Aug 2016 at 17:37 Onno Bokhove <O.Bokhove(a)leeds.ac.uk> >>>>> wrote: >>>>> >>>>> It is possible to merge the branch for the mixed system >>>>> relatively soon? >>>>> Anna's results are relevant for the grant application Colin and I aim >>>>> to submit in Sept. on further wave-ship modelling. >>>>> Sorry for asking! >>>>> ------------------------------ >>>>> *From:* firedrake-bounces(a)imperial.ac.uk <firedrake- >>>>> bounces(a)imperial.ac.uk> on behalf of David Ham < >>>>> David.Ham(a)imperial.ac.uk> >>>>> *Sent:* Friday, August 5, 2016 11:37:23 AM >>>>> *To:* firedrake >>>>> *Subject:* Re: [firedrake] Mixed system >>>>> >>>>> >>>>> >>>>> On Fri, 5 Aug 2016 at 11:31 Anna Kalogirou <a.kalogirou(a)leeds.ac.uk> >>>>> wrote: >>>>> >>>>> Hi Kramer, >>>>> >>>>> Thank you for the reply. >>>>> >>>>> The only problem I see with the solution you suggest, is that both I >>>>> and >>>>> lambda will be trial functions, so both terms should be on the bilinear >>>>> part of the variational form. >>>>> >>>>> >>>>> Not quite. I and lambda together form a Function in a mixed finite >>>>> element space. You solve for both of them together, just like pressure and >>>>> velocity in a fluid simulation. >>>>> >>>>> >>>>> >>>>> I could, of course, define a residual and solve everything using >>>>> nonlinear solvers, but I don't know how efficient that would be. >>>>> >>>>> >>>>> I think you are confused. None of the relevant issues change in the >>>>> nonlinear case. After all, the nonlinear solver will differentiate the >>>>> residual to create a bilinear form anyway. >>>>> >>>>> Incidentally, there is no particular reason to believe that writing >>>>> the problem in nonlinear form would be inefficient. >>>>> >>>>> Regards, >>>>> >>>>> David >>>>> >>>>> >>>>> >>>>> Best, >>>>> >>>>> Anna. >>>>> >>>>> Dr Anna Kalogirou >>>>> Research Fellow >>>>> School of Mathematics >>>>> University of Leeds >>>>> >>>>> http://www1.maths.leeds.ac.uk/~matak/ >>>>> >>>>> On 04/08/16 17:29, Stephan Kramer wrote: >>>>> > On 27/07/16 11:01, Anna Kalogirou wrote: >>>>> >> Dear all, >>>>> >> >>>>> >> I am trying to solve the system found in the attached pdf as a mixed >>>>> >> system. I want to solve the same >>>>> >> problem as the one found here >>>>> >> <https://bitbucket.org/annakalog/buoy2d/src/ >>>>> 14334d3c20b9f10ed7c1246cde9e3cb60b1c75e4/Inequality% >>>>> 20constraint/?at=master>, >>>>> >> >>>>> >> but with the use of Schur complements and not linear algebra. >>>>> >> >>>>> >> In a previous discussion about this problem, Lawrence mentioned >>>>> that the >>>>> >> test function for integral(lambda*Theta(x-Lp)dx) needs to be >>>>> considered >>>>> >> as coming from the real space of constant functions. Could you >>>>> please >>>>> >> elaborate on this? >>>>> > >>>>> > My guess at what Lawrence means is (it's how I would probably >>>>> approach >>>>> > it), is that you add an additional equation that says >>>>> > >>>>> > (1) I = integral(lambda*Theta(x-Lp)dx) >>>>> > >>>>> > and then substitute I in your third equation. Equation (1) can >>>>> > be turned into a finite element weak formulation by considering the >>>>> > function space of functions that are constant over the entire domain >>>>> > This function space is just 1-dimensional and therefore equivalent to >>>>> > the space of reals R and hence we denote this function >>>>> >>>>> -- http://www.imperial.ac.uk/people/colin.cotter www.cambridge.org/9781107663916

3 4

Re: [firedrake] Mixed system
by David Ham 05 Oct '16

05 Oct '16

Hi Anna, Are you certain that the test you are running is the one you pointed to on the web. Line 89 from your error message does not appear in the version on bitbucket, and the last commit was 4 days ago. Regards, David On Fri, 16 Sep 2016 at 17:03 Anna Kalogirou <A.Kalogirou(a)leeds.ac.uk> wrote: > When I run > > firedrake-zenodo -m "Versions I am using now” > > I get > > Failed to retrieve git hash for optional component 'SLOPE', continuing > without it > > The file firedrake.json contains the following (up to ufl it’s the same > as yours): > > {"fiat": "1bc440ca34d4b127897cf5a36d63bf8a41a3d55c", "COFFEE": > "8e03328d659fd8d8871bf5b7bc3f08b769f0ab9d", "ufl": > "2fc6d00f7e4b5e0639c6d01f27f9b2bbc0d4696e", "firedrake": > "71cf933f05508427287dc6464e15d2d7adf1b179", "tsfc": > "f867938c58f475f1cd332421e2bf14cffa2f5aef", "petsc4py": > "0b590afa9d3c42e0b4f863061b6e082d6f0227f6", "petsc": > "21143a8c817f512a23f8e18da648891eb40f5217", "PyOP2": > "d9ea13e1e17fcb1e05adcf1db31bd8190974520d", "message": "Versions I am using > now"} > > > On 16 Sep 2016, at 16:49, David Ham <David.Ham(a)imperial.ac.uk> wrote: > > OK. Let's do some sanity checking. > > 1. Please confirm that the current firedrake checkout is in developer mode > and has always been in developer mode (I can see that it is in developer > mode from the error message but I want to make sure there cannot be stale > files around from a non-developer mode past). > > 2. Run firedrake-clean and check that the error still occurs. > > 2. Please run the following command: > > firedrake-zenodo -m "Versions I am using now" > > This will create a file firedrake.json . Please check that its contents > exactly match what is below: > > {"fiat": "1bc440ca34d4b127897cf5a36d63bf8a41a3d55c", "COFFEE": > "8e03328d659fd8d8871bf5b7bc3f08b769f0ab9d", "ufl": > "2fc6d00f7e4b5e0639c6d01f27f9b2bbc0d4696e", "firedrake": > "71cf933f05508427287dc6464e15d2d7adf1b179", "tsfc": > "f867938c58f475f1cd332421e2bf14cffa2f5aef", "petsc4py": > "0b590afa9d3c42e0b4f863061b6e082d6f0227f6", "petsc": > "21143a8c817f512a23f8e18da648891eb40f5217", "PyOP2": > "d9ea13e1e17fcb1e05adcf1db31bd8190974520d", "message": "Versions I am using > now"} > > > On Fri, 16 Sep 2016 at 16:28 Anna Kalogirou <A.Kalogirou(a)leeds.ac.uk> > wrote: > >> That is the same version I have. >> >> >> On 16 Sep 2016, at 16:23, David Ham <David.Ham(a)imperial.ac.uk> wrote: >> >> Hi Anna, >> >> Please confirm the git revision of your PyOP2. You can retrieve this by >> typing the following command in the PyOP2 directory: >> >> git rev-parse HEAD >> >> You should have: >> d9ea13e1e17fcb1e05adcf1db31bd8190974520d >> >> Regards, >> >> David >> >> On Fri, 16 Sep 2016 at 15:59 Anna Kalogirou <A.Kalogirou(a)leeds.ac.uk> >> wrote: >> >>> Hi David, >>> >>> Thank you for the reply. Indeed, the problem was with the virtualenv; I >>> had to deactivate and activate it again. Unfortunately I still get errors, >>> which are probably due to the firedrake I have installed (since you said >>> that it works on your computer): >>> >>> File "buoy-swe.py", line 89, in <module> >>> F_solver = solver_F(phi0_5, eta1, mu0_5, I, w, phi0, eta0, Z0, W0, >>> step_b, etaR, phi_t, eta_t, mu_t, I_t, v1, v2, v3, v4, dt, Hb, H0, L, >>> dR_dt, g, rho, Mass, solvers_print); >>> File "/Users/matak/Documents/Simulations/Firedrake/Ship/Modules/Mixed >>> system/solvers.py", line 41, in solver_F >>> F_solver = LinearVariationalSolver(F_problem, >>> solver_parameters=solvers_print) >>> File >>> "/Users/matak/firedrake/src/firedrake/firedrake/variational_solver.py", >>> line 264, in __init__ >>> super(LinearVariationalSolver, self).__init__(*args, **kwargs) >>> File >>> "/Users/matak/firedrake/src/firedrake/firedrake/variational_solver.py", >>> line 127, in __init__ >>> ctx = solving_utils._SNESContext(problem) >>> File >>> "/Users/matak/firedrake/src/firedrake/firedrake/solving_utils.py", line >>> 107, in __init__ >>> for problem in problems) >>> File >>> "/Users/matak/firedrake/src/firedrake/firedrake/solving_utils.py", line >>> 107, in <genexpr> >>> for problem in problems) >>> File "/Users/matak/firedrake/src/firedrake/firedrake/assemble.py", >>> line 67, in assemble >>> inverse=inverse, nest=nest) >>> File "<decorator-gen-294>", line 2, in _assemble >>> File "/Users/matak/firedrake/src/firedrake/firedrake/utils.py", line >>> 62, in wrapper >>> return f(*args, **kwargs) >>> File "/Users/matak/firedrake/src/firedrake/firedrake/assemble.py", >>> line 179, in _assemble >>> nest=nest) >>> File "/Users/matak/firedrake/src/PyOP2/pyop2/backends.py", line 118, >>> in __call__ >>> return t(*args, **kwargs) >>> File "/Users/matak/firedrake/src/PyOP2/pyop2/caching.py", line 122, in >>> __new__ >>> args, kwargs = cls._process_args(*args, **kwargs) >>> File "<decorator-gen-258>", line 2, in _process_args >>> File "/Users/matak/firedrake/src/PyOP2/pyop2/utils.py", line 130, in >>> wrapper >>> return f(*args, **kwargs) >>> File "/Users/matak/firedrake/src/PyOP2/pyop2/base.py", line 3591, in >>> _process_args >>> if not (pair[0].toset == dsets[0].set and >>> File "/Users/matak/firedrake/src/PyOP2/pyop2/utils.py", line 64, in >>> __get__ >>> obj.__dict__[self.__name__] = result = self.fget(obj) >>> File "/Users/matak/firedrake/src/PyOP2/pyop2/base.py", line 3372, in >>> toset >>> m.toset for m in self._maps)) >>> File "/Users/matak/firedrake/src/PyOP2/pyop2/caching.py", line 156, in >>> __new__ >>> obj = make_obj() >>> File "/Users/matak/firedrake/src/PyOP2/pyop2/caching.py", line 137, in >>> make_obj >>> obj.__init__(*args, **kwargs) >>> File "/Users/matak/firedrake/src/PyOP2/pyop2/base.py", line 977, in >>> __init__ >>> "All components of a MixedSet must have the same number of layers." >>> AssertionError: All components of a MixedSet must have the same number >>> of layers. >>> >>> >>> >>> On 16 Sep 2016, at 15:47, David Ham <David.Ham(a)imperial.ac.uk> wrote: >>> >>> Hi Anna, >>> >>> Your code works for me (in the sense that it appears to run without >>> error. I have no idea whether the answer is right ;) >>> >>> Please also note that Miklos is correct: your error is definitely from a >>> global firedrake install. You can tell because the paths start /usr/local. >>> Your previous error messages are from a local install. You can tell that >>> because the paths start /Users/matak . Are you sure your virtualenv is >>> active? >>> >>> Regards, >>> >>> David >>> >>> >>> >>> On Fri, 16 Sep 2016 at 15:16 Anna Kalogirou <A.Kalogirou(a)leeds.ac.uk> >>> wrote: >>> >>>> During the past few days I have been working on my laptop (not the >>>> desktop at Leeds) and I don’t believe there are multiple installations of >>>> firedrake on this machine. >>>> >>>> Can someone please try to run my code >>>> <https://bitbucket.org/annakalog/buoy2d/src/673fef3466b06c98ac11f5a120ff4211…> so >>>> that we know for sure if the problem is in my code or related with the >>>> installation? Thanks. >>>> >>>> >>>> On 16 Sep 2016, at 14:54, Colin Cotter <colin.cotter(a)imperial.ac.uk> >>>> wrote: >>>> >>>> Ah yes, that tallies with Daniel Ruprecht's reports that Firedrake is >>>> part of the standard installation at Leeds! >>>> >>>> On 16 September 2016 at 14:53, Homolya, Miklós < >>>> m.homolya14(a)imperial.ac.uk> wrote: >>>> >>>> It seems that the Firedrake you're importing is a system wide >>>>> installation. No wonder that updating your own installation doesn't fix it. >>>>> >>>>> Get Outlook for Android <https://aka.ms/ghei36> >>>>> >>>> >>>>> >>>>> >>>>> On Thu, Sep 15, 2016 at 8:33 PM +0100, "Anna Kalogirou" < >>>>> A.Kalogirou(a)leeds.ac.uk> wrote: >>>>> >>>>> Ok thanks, I get the error below: >>>>> >>>>> Traceback (most recent call last): >>>>> >>>> File "buoy-swe.py", line 6, in <module> >>>>> from firedrake import * >>>>> File "/usr/local/lib/python2.7/site-packages/firedrake/__init__.py", >>>>> line 30, in <module> >>>>> from assemble import * >>>>> File "/usr/local/lib/python2.7/site-packages/firedrake/assemble.py", >>>>> line 8, in <module> >>>>> import assembly_cache >>>>> File >>>>> "/usr/local/lib/python2.7/site-packages/firedrake/assembly_cache.py", line >>>>> 40, in <module> >>>>> import function >>>>> File "/usr/local/lib/python2.7/site-packages/firedrake/function.py", >>>>> line 9, in <module> >>>>> import assemble_expressions >>>>> File >>>>> "/usr/local/lib/python2.7/site-packages/firedrake/assemble_expressions.py", >>>>> line 15, in <module> >>>>> import functionspace >>>>> File >>>>> "/usr/local/lib/python2.7/site-packages/firedrake/functionspace.py", line >>>>> 10, in <module> >>>>> import dmplex >>>>> ImportError: >>>>> dlopen(/usr/local/lib/python2.7/site-packages/firedrake/dmplex.so, 2): >>>>> Library not loaded: >>>>> /usr/local/lib/python2.7/site-packages/petsc/lib/libpetsc.3.05.2.dylib >>>>> Referenced from: >>>>> /usr/local/lib/python2.7/site-packages/firedrake/dmplex.so >>>>> Reason: image not found >>>>> >>>> On 15 Sep 2016, at 18:21, Andrew McRae <A.T.T.McRae(a)bath.ac.uk> wrote: >>>>> >>>>> To give credit where it's due, >>>>> >>>>> s/Andrew's patch/Miklós's patch/ >>>>> >>>>> On 15 September 2016 at 18:16, David Ham <David.Ham(a)imperial.ac.uk> >>>>> wrote: >>>>> >>>>> I suspect Andrew is correct. However the deeper issue is that you >>>>>> edited the installed file, which is never a safe thing to do. One must >>>>>> always edit the file in the source directory and then (if necessary) >>>>>> re-install. >>>>>> >>>>>> Anyway, I have committed Andrew's patch to the repository, so >>>>>> hopefully the problem will go away if you do a git update. >>>>>> >>>>>> Regards, >>>>>> >>>>>> David >>>>>> >>>>>> On Wed, 14 Sep 2016 at 18:39 Andrew McRae <A.T.T.McRae(a)bath.ac.uk> >>>>>> wrote: >>>>>> >>>>>>> Probably related to installing Firedrake in development mode or not. >>>>>>> >>>>>>> I can't help any further here; I don't know if someone else can. >>>>>>> >>>>>>> On 14 September 2016 at 16:13, Anna Kalogirou < >>>>>>> A.Kalogirou(a)leeds.ac.uk> wrote: >>>>>>> >>>>>>>> Ok, so I replace line 3366 in file “ >>>>>>>> /Users/matak/firedrake/lib/python2.7/site-packages/pyop2/base.py” >>>>>>>> with the suggested patch, but I now get the error below. >>>>>>>> >>>>>>>> Note that the error I get comes from file >>>>>>>> firedrake/lib/python2.7/site-packages/pyop2/base.py >>>>>>>> while in the reported issue the error comes from >>>>>>>> firedrake/src/PyOP2/pyop2/base.py >>>>>>>> so the two might not be related after all. >>>>>>>> >>>>>>>> Thanks, >>>>>>>> Anna. >>>>>>>> >>>>>>> >>>>>>>> >>>>>>>> Traceback (most recent call last): >>>>>>>> File "buoy-swe.py", line 89, in <module> >>>>>>>> F_solver = solver_F(phi0_5, eta1, mu0_5, I, w, phi0, eta0, Z0, >>>>>>>> W0, step_b, etaR, phi_t, eta_t, mu_t, I_t, v1, v2, v3, v4, dt, Hb, H0, L, >>>>>>>> dR_dt, g, rho, Mass, solvers_print); >>>>>>>> File >>>>>>>> "/Users/matak/Documents/Simulations/Firedrake/Ship/Modules/Mixed >>>>>>>> system/solvers.py", line 41, in solver_F >>>>>>>> F_solver = LinearVariationalSolver(F_problem, >>>>>>>> solver_parameters=solvers_print) >>>>>>>> File >>>>>>>> "/Users/matak/firedrake/lib/python2.7/site-packages/firedrake/variational_solver.py", >>>>>>>> line 264, in __init__ >>>>>>>> super(LinearVariationalSolver, self).__init__(*args, **kwargs) >>>>>>>> File >>>>>>>> "/Users/matak/firedrake/lib/python2.7/site-packages/firedrake/variational_solver.py", >>>>>>>> line 127, in __init__ >>>>>>>> ctx = solving_utils._SNESContext(problem) >>>>>>>> File >>>>>>>> "/Users/matak/firedrake/lib/python2.7/site-packages/firedrake/solving_utils.py", >>>>>>>> line 107, in __init__ >>>>>>>> for problem in problems) >>>>>>>> File >>>>>>>> "/Users/matak/firedrake/lib/python2.7/site-packages/firedrake/solving_utils.py", >>>>>>>> line 107, in <genexpr> >>>>>>>> for problem in problems) >>>>>>>> File >>>>>>>> "/Users/matak/firedrake/lib/python2.7/site-packages/firedrake/assemble.py", >>>>>>>> line 67, in assemble >>>>>>>> inverse=inverse, nest=nest) >>>>>>>> File "<decorator-gen-294>", line 2, in _assemble >>>>>>>> File >>>>>>>> "/Users/matak/firedrake/lib/python2.7/site-packages/firedrake/utils.py", >>>>>>>> line 62, in wrapper >>>>>>>> return f(*args, **kwargs) >>>>>>>> File >>>>>>>> "/Users/matak/firedrake/lib/python2.7/site-packages/firedrake/assemble.py", >>>>>>>> line 179, in _assemble >>>>>>>> nest=nest) >>>>>>>> File >>>>>>>> "/Users/matak/firedrake/lib/python2.7/site-packages/pyop2/backends.py", >>>>>>>> line 118, in __call__ >>>>>>>> return t(*args, **kwargs) >>>>>>>> File >>>>>>>> "/Users/matak/firedrake/lib/python2.7/site-packages/pyop2/caching.py", line >>>>>>>> 122, in __new__ >>>>>>>> args, kwargs = cls._process_args(*args, **kwargs) >>>>>>>> File "<decorator-gen-258>", line 2, in _process_args >>>>>>>> File >>>>>>>> "/Users/matak/firedrake/lib/python2.7/site-packages/pyop2/utils.py", line >>>>>>>> 130, in wrapper >>>>>>>> return f(*args, **kwargs) >>>>>>>> File >>>>>>>> "/Users/matak/firedrake/lib/python2.7/site-packages/pyop2/base.py", line >>>>>>>> 3592, in _process_args >>>>>>>> if not (pair[0].toset == dsets[0].set and >>>>>>>> File >>>>>>>> "/Users/matak/firedrake/lib/python2.7/site-packages/pyop2/utils.py", line >>>>>>>> 64, in __get__ >>>>>>>> obj.__dict__[self.__name__] = result = self.fget(obj) >>>>>>>> File >>>>>>>> "/Users/matak/firedrake/lib/python2.7/site-packages/pyop2/base.py", line >>>>>>>> 3373, in toset >>>>>>>> m.toset for m in self._maps)) >>>>>>>> File >>>>>>>> "/Users/matak/firedrake/lib/python2.7/site-packages/pyop2/caching.py", line >>>>>>>> 156, in __new__ >>>>>>>> obj = make_obj() >>>>>>>> File >>>>>>>> "/Users/matak/firedrake/lib/python2.7/site-packages/pyop2/caching.py", line >>>>>>>> 137, in make_obj >>>>>>>> obj.__init__(*args, **kwargs) >>>>>>>> File >>>>>>>> "/Users/matak/firedrake/lib/python2.7/site-packages/pyop2/base.py", line >>>>>>>> 977, in __init__ >>>>>>>> "All components of a MixedSet must have the same number of >>>>>>>> layers." >>>>>>>> AssertionError: All components of a MixedSet must have the same >>>>>>>> number of layers. >>>>>>>> >>>>>>> >>>>>>>> ********************************************************** >>>>>>>> >>>>>>> >>>>>>>> Dr Anna Kalogirou >>>>>>>> Research Fellow >>>>>>>> >>>>>>> Department of Applied Mathematics >>>>>>>> University of Leeds >>>>>>>> >>>>>>>> Web: http://www1.maths.leeds.ac.uk/~matak/ >>>>>>>> >>>>>>> Youtube: https://www.youtube.com/channel/UCUZAYHtVoiMqepQflisp66g >>>>>>>> >>>>>>>> ********************************************************** >>>>>>>> Leverhulme Trust Early Career Fellow 2017-2019 >>>>>>>> >>>>>>>> ********************************************************** >>>>>>>> >>>>>>>> >>>>>>>> On 14 Sep 2016, at 16:00, Andrew McRae <A.T.T.McRae(a)bath.ac.uk> >>>>>>>> wrote: >>>>>>>> >>>>>>>> Not as a branch, but you can change PyOP2 locally using the code at >>>>>>>> the bottom of that post, "Suggested patch" >>>>>>>> >>>>>>>> On 14 September 2016 at 15:54, Anna Kalogirou < >>>>>>>> A.Kalogirou(a)leeds.ac.uk> wrote: >>>>>>>> >>>>>>>> Hi Andrew, >>>>>>>>> >>>>>>>>> Yes indeed, the two might be related. Was there a fix for the >>>>>>>>> issue reported on github? >>>>>>>>> >>>>>>>>> Best, Anna. >>>>>>>>> >>>>>>>>> >>>>>>>>> ********************************************************** >>>>>>>>> >>>>>>>> >>>>>>>>> Dr Anna Kalogirou >>>>>>>>> Research Fellow >>>>>>>>> >>>>>>>> Department of Applied Mathematics >>>>>>>>> University of Leeds >>>>>>>>> >>>>>>>>> Web: http://www1.maths.leeds.ac.uk/~matak/ >>>>>>>>> >>>>>>>> Youtube: https://www.youtube.com/channel/UCUZAYHtVoiMqepQflisp66g >>>>>>>>> >>>>>>>>> ********************************************************** >>>>>>>>> Leverhulme Trust Early Career Fellow 2017-2019 >>>>>>>>> >>>>>>>>> ********************************************************** >>>>>>>>> >>>>>>>>> >>>>>>>>> On 13 Sep 2016, at 16:34, Andrew McRae <A.T.T.McRae(a)bath.ac.uk> >>>>>>>>> wrote: >>>>>>>>> >>>>>>>>> Probably related to >>>>>>>>> https://github.com/firedrakeproject/firedrake/issues/862 (third >>>>>>>>> post, not first). You might be able to adapt the fix; I haven't looked in >>>>>>>>> detail. >>>>>>>>> >>>>>>>>> On 12 September 2016 at 16:34, Anna Kalogirou < >>>>>>>>> A.Kalogirou(a)leeds.ac.uk> wrote: >>>>>>>>> >>>>>>>>> Dear all, >>>>>>>>>> >>>>>>>>>> Regarding the problem described in my email sent of Friday, I >>>>>>>>>> fixed a minor typo and now get the error below. Any ideas what might be >>>>>>>>>> causing it? Thanks! >>>>>>>>>> >>>>>>>>>> Best, Anna. >>>>>>>>>> >>>>>>>>>> Traceback (most recent call last): >>>>>>>>>> File "buoy-swe.py", line 89, in <module> >>>>>>>>>> F_solver = solver_F(phi0_5, eta1, mu0_5, I, w, phi0, eta0, >>>>>>>>>> Z0, W0, step_b, etaR, phi_t, eta_t, mu_t, I_t, v1, v2, v3, v4, dt, Hb, H0, >>>>>>>>>> L, dR_dt, g, rho, Mass, solvers_print); >>>>>>>>>> File >>>>>>>>>> "/Users/matak/Documents/Simulations/Firedrake/Ship/Modules/Mixed >>>>>>>>>> system/solvers.py", line 41, in solver_F >>>>>>>>>> F_solver = LinearVariationalSolver(F_problem, >>>>>>>>>> solver_parameters=solvers_print) >>>>>>>>>> File >>>>>>>>>> "/Users/matak/firedrake/lib/python2.7/site-packages/firedrake/variational_solver.py", >>>>>>>>>> line 264, in __init__ >>>>>>>>>> super(LinearVariationalSolver, self).__init__(*args, **kwargs) >>>>>>>>>> File >>>>>>>>>> "/Users/matak/firedrake/lib/python2.7/site-packages/firedrake/variational_solver.py", >>>>>>>>>> line 127, in __init__ >>>>>>>>>> ctx = solving_utils._SNESContext(problem) >>>>>>>>>> File >>>>>>>>>> "/Users/matak/firedrake/lib/python2.7/site-packages/firedrake/solving_utils.py", >>>>>>>>>> line 107, in __init__ >>>>>>>>>> for problem in problems) >>>>>>>>>> File >>>>>>>>>> "/Users/matak/firedrake/lib/python2.7/site-packages/firedrake/solving_utils.py", >>>>>>>>>> line 107, in <genexpr> >>>>>>>>>> for problem in problems) >>>>>>>>>> File >>>>>>>>>> "/Users/matak/firedrake/lib/python2.7/site-packages/firedrake/assemble.py", >>>>>>>>>> line 67, in assemble >>>>>>>>>> inverse=inverse, nest=nest) >>>>>>>>>> File "<decorator-gen-294>", line 2, in _assemble >>>>>>>>>> File >>>>>>>>>> "/Users/matak/firedrake/lib/python2.7/site-packages/firedrake/utils.py", >>>>>>>>>> line 62, in wrapper >>>>>>>>>> return f(*args, **kwargs) >>>>>>>>>> File >>>>>>>>>> "/Users/matak/firedrake/lib/python2.7/site-packages/firedrake/assemble.py", >>>>>>>>>> line 163, in _assemble >>>>>>>>>> map_pairs.append((op2.DecoratedMap(test.cell_node_map(), >>>>>>>>>> cell_domains), >>>>>>>>>> File >>>>>>>>>> "/Users/matak/firedrake/lib/python2.7/site-packages/firedrake/functionspaceimpl.py", >>>>>>>>>> line 646, in cell_node_map >>>>>>>>>> for i, s in enumerate(self._spaces)) >>>>>>>>>> File >>>>>>>>>> "/Users/matak/firedrake/lib/python2.7/site-packages/pyop2/backends.py", >>>>>>>>>> line 118, in __call__ >>>>>>>>>> return t(*args, **kwargs) >>>>>>>>>> File >>>>>>>>>> "/Users/matak/firedrake/lib/python2.7/site-packages/pyop2/caching.py", line >>>>>>>>>> 156, in __new__ >>>>>>>>>> obj = make_obj() >>>>>>>>>> File >>>>>>>>>> "/Users/matak/firedrake/lib/python2.7/site-packages/pyop2/caching.py", line >>>>>>>>>> 137, in make_obj >>>>>>>>>> obj.__init__(*args, **kwargs) >>>>>>>>>> File >>>>>>>>>> "/Users/matak/firedrake/lib/python2.7/site-packages/pyop2/base.py", line >>>>>>>>>> 3342, in __init__ >>>>>>>>>> if not all(m is None or m.iterset == self.iterset for m in >>>>>>>>>> self._maps): >>>>>>>>>> File >>>>>>>>>> "/Users/matak/firedrake/lib/python2.7/site-packages/pyop2/base.py", line >>>>>>>>>> 3342, in <genexpr> >>>>>>>>>> if not all(m is None or m.iterset == self.iterset for m in >>>>>>>>>> self._maps): >>>>>>>>>> File >>>>>>>>>> "/Users/matak/firedrake/lib/python2.7/site-packages/pyop2/utils.py", line >>>>>>>>>> 64, in __get__ >>>>>>>>>> obj.__dict__[self.__name__] = result = self.fget(obj) >>>>>>>>>> File >>>>>>>>>> "/Users/matak/firedrake/lib/python2.7/site-packages/pyop2/base.py", line >>>>>>>>>> 3366, in iterset >>>>>>>>>> return reduce(lambda a, b: a if a is None else a.iterset or b >>>>>>>>>> if b is None else b.iterset, self._maps) >>>>>>>>>> File >>>>>>>>>> "/Users/matak/firedrake/lib/python2.7/site-packages/pyop2/base.py", line >>>>>>>>>> 3366, in <lambda> >>>>>>>>>> return reduce(lambda a, b: a if a is None else a.iterset or b >>>>>>>>>> if b is None else b.iterset, self._maps) >>>>>>>>>> File >>>>>>>>>> "/Users/matak/firedrake/lib/python2.7/site-packages/pyop2/base.py", line >>>>>>>>>> 799, in __getattr__ >>>>>>>>>> return getattr(self._parent, name) >>>>>>>>>> AttributeError: 'Set' object has no attribute ‘iterset' >>>>>>>>>> >>>>>>>>>> >>>>>>>>>> ********************************************************** >>>>>>>>>> >>>>>>>>> >>>>>>>>>> Dr Anna Kalogirou >>>>>>>>>> Research Fellow >>>>>>>>>> >>>>>>>>> Department of Applied Mathematics >>>>>>>>>> University of Leeds >>>>>>>>>> >>>>>>>>>> Web: http://www1.maths.leeds.ac.uk/~matak/ >>>>>>>>>> >>>>>>>>> Youtube: https://www.youtube.com/channel/UCUZAYHtVoiMqepQflisp66g >>>>>>>>>> >>>>>>>>>> ********************************************************** >>>>>>>>>> Leverhulme Trust Early Career Fellow 2017-2019 >>>>>>>>>> >>>>>>>>>> ********************************************************** >>>>>>>>>> >>>>>>>>>> >>>>>>>>>> On 9 Sep 2016, at 17:37, Anna Kalogirou <a.kalogirou(a)leeds.ac.uk> >>>>>>>>>> wrote: >>>>>>>>>> >>>>>>>>>> Hi David, >>>>>>>>>> >>>>>>>>>> I have tried to implement my problem (attached) using the real >>>>>>>>>> function spaces but I get the error below. The code can be found >>>>>>>>>> here >>>>>>>>>> <https://bitbucket.org/annakalog/buoy2d/src/d6a805b6a8d0/Mixed%20system/?at=…> >>>>>>>>>> . >>>>>>>>>> >>>>>>>>>> I define the real space by V_R = FunctionSpace(mesh, "R", 0). Is >>>>>>>>>> that correct? >>>>>>>>>> Also, is it a problem that the linear part of the variational >>>>>>>>>> form does not depend on one of the test functions from the mixed function >>>>>>>>>> space (only depends on v1, v2, v3 and not v4)? >>>>>>>>>> >>>>>>>>>> Best, >>>>>>>>>> >>>>>>>>>> Anna. >>>>>>>>>> >>>>>>>>>> >>>>>>>>>> >>>>>>>>>> Traceback (most recent call last): >>>>>>>>>> File "buoy-swe.py", line 89, in <module> >>>>>>>>>> F_solver = solver_F(phi0_5, eta1, mu0_5, I, w, phi0, eta0, >>>>>>>>>> Z0, W0, step_b, etaR, phi_t, eta_t, mu_t, I_t, v1, v2, v3, v4, dt, Hb, H0, >>>>>>>>>> L, dR_dt, g, rho, Mass, solvers_print); >>>>>>>>>> File >>>>>>>>>> "/Users/matak/Documents/Simulations/Firedrake/Ship/Modules/Mixed >>>>>>>>>> system/solvers.py", line 41, in solver_F >>>>>>>>>> F_solver = LinearVariationalSolver(F_problem, >>>>>>>>>> solver_parameters=solvers_print) >>>>>>>>>> File >>>>>>>>>> "/Users/matak/firedrake/lib/python2.7/site-packages/firedrake/variational_solver.py", >>>>>>>>>> line 264, in __init__ >>>>>>>>>> super(LinearVariationalSolver, self).__init__(*args, **kwargs) >>>>>>>>>> File >>>>>>>>>> "/Users/matak/firedrake/lib/python2.7/site-packages/firedrake/variational_solver.py", >>>>>>>>>> line 127, in __init__ >>>>>>>>>> ctx = solving_utils._SNESContext(problem) >>>>>>>>>> File >>>>>>>>>> "/Users/matak/firedrake/lib/python2.7/site-packages/firedrake/solving_utils.py", >>>>>>>>>> line 107, in __init__ >>>>>>>>>> for problem in problems) >>>>>>>>>> File >>>>>>>>>> "/Users/matak/firedrake/lib/python2.7/site-packages/firedrake/solving_utils.py", >>>>>>>>>> line 107, in <genexpr> >>>>>>>>>> for problem in problems) >>>>>>>>>> File >>>>>>>>>> "/Users/matak/firedrake/lib/python2.7/site-packages/firedrake/assemble.py", >>>>>>>>>> line 67, in assemble >>>>>>>>>> inverse=inverse, nest=nest) >>>>>>>>>> File "<decorator-gen-294>", line 2, in _assemble >>>>>>>>>> File >>>>>>>>>> "/Users/matak/firedrake/lib/python2.7/site-packages/firedrake/utils.py", >>>>>>>>>> line 62, in wrapper >>>>>>>>>> return f(*args, **kwargs) >>>>>>>>>> File >>>>>>>>>> "/Users/matak/firedrake/lib/python2.7/site-packages/firedrake/assemble.py", >>>>>>>>>> line 101, in _assemble >>>>>>>>>> inverse=inverse) >>>>>>>>>> File >>>>>>>>>> "/Users/matak/firedrake/lib/python2.7/site-packages/firedrake/tsfc_interface.py", >>>>>>>>>> line 307, in compile_form >>>>>>>>>> number_map).kernels >>>>>>>>>> File >>>>>>>>>> "/Users/matak/firedrake/lib/python2.7/site-packages/pyop2/caching.py", line >>>>>>>>>> 198, in __new__ >>>>>>>>>> obj = make_obj() >>>>>>>>>> File >>>>>>>>>> "/Users/matak/firedrake/lib/python2.7/site-packages/pyop2/caching.py", line >>>>>>>>>> 188, in make_obj >>>>>>>>>> obj.__init__(*args, **kwargs) >>>>>>>>>> File >>>>>>>>>> "/Users/matak/firedrake/lib/python2.7/site-packages/firedrake/tsfc_interface.py", >>>>>>>>>> line 212, in __init__ >>>>>>>>>> tree = tsfc_compile_form(form, prefix=name, >>>>>>>>>> parameters=parameters) >>>>>>>>>> File >>>>>>>>>> "/Users/matak/firedrake/lib/python2.7/site-packages/tsfc/driver.py", line >>>>>>>>>> 47, in compile_form >>>>>>>>>> do_estimate_degrees=True) >>>>>>>>>> File >>>>>>>>>> "/Users/matak/firedrake/lib/python2.7/site-packages/ufl/algorithms/compute_form_data.py", >>>>>>>>>> line 386, in compute_form_data >>>>>>>>>> check_form_arity(preprocessed_form, >>>>>>>>>> self.original_form.arguments()) # Currently testing how fast this is >>>>>>>>>> File >>>>>>>>>> "/Users/matak/firedrake/lib/python2.7/site-packages/ufl/algorithms/check_arities.py", >>>>>>>>>> line 152, in check_form_arity >>>>>>>>>> check_integrand_arity(itg.integrand(), arguments) >>>>>>>>>> File >>>>>>>>>> "/Users/matak/firedrake/lib/python2.7/site-packages/ufl/algorithms/check_arities.py", >>>>>>>>>> line 145, in check_integrand_arity >>>>>>>>>> args = map_expr_dag(rules, expr, compress=False) >>>>>>>>>> File >>>>>>>>>> "/Users/matak/firedrake/lib/python2.7/site-packages/ufl/corealg/map_dag.py", >>>>>>>>>> line 37, in map_expr_dag >>>>>>>>>> result, = map_expr_dags(function, [expression], >>>>>>>>>> compress=compress) >>>>>>>>>> File >>>>>>>>>> "/Users/matak/firedrake/lib/python2.7/site-packages/ufl/corealg/map_dag.py", >>>>>>>>>> line 86, in map_expr_dags >>>>>>>>>> r = handlers[v._ufl_typecode_](v, *[vcache[u] for u in >>>>>>>>>> v.ufl_operands]) >>>>>>>>>> File >>>>>>>>>> "/Users/matak/firedrake/lib/python2.7/site-packages/ufl/algorithms/check_arities.py", >>>>>>>>>> line 42, in sum >>>>>>>>>> raise ArityMismatch("Adding expressions with non-matching >>>>>>>>>> form arguments {0} vs {1}.".format(a, b)) >>>>>>>>>> ufl.algorithms.check_arities.ArityMismatch: Adding expressions >>>>>>>>>> with non-matching form arguments () vs >>>>>>>>>> (Argument(WithGeometry(IndexedProxyFunctionSpace(<firedrake.mesh.ExtrudedMeshTopology >>>>>>>>>> object at 0x113e6b890>, TensorProductElement(FiniteElement('Lagrange', >>>>>>>>>> interval, 1), FiniteElement('Lagrange', interval, 1), >>>>>>>>>> cell=TensorProductCell(interval, interval)), name=None, index=0, >>>>>>>>>> component=None), >>>>>>>>>> Mesh(VectorElement(TensorProductElement(FiniteElement('Lagrange', interval, >>>>>>>>>> 1), FiniteElement('Lagrange', interval, 1), >>>>>>>>>> cell=TensorProductCell(interval, interval)), dim=2), 4)), 1, None),). >>>>>>>>>> <coupled_system.pdf> >>>>>>>>>> >>>>>>>>>> On 18 Aug 2016, at 16:03, David Ham <david.ham(a)imperial.ac.uk> >>>>>>>>>> wrote: >>>>>>>>>> >>>>>>>>>> Hi Anna, >>>>>>>>>> >>>>>>>>>> Real function spaces now pass some simple tests, at least on my >>>>>>>>>> machine. You need the real-function-space branch of firedrake, and the >>>>>>>>>> globals_matrices branch of PyOP2. It's not going to get merged before I >>>>>>>>>> disappear for a week - not least because I need to update documentation. >>>>>>>>>> However you should be able to test it out on the branches. >>>>>>>>>> >>>>>>>>>> Regards, >>>>>>>>>> >>>>>>>>>> David >>>>>>>>>> >>>>>>>>>> On Wed, 10 Aug 2016 at 16:58 David Ham <David.Ham(a)imperial.ac.uk> >>>>>>>>>> wrote: >>>>>>>>>> >>>>>>>>>>> Unfortunately, It isn't as simple as saying "no mpi". MPI >>>>>>>>>>> pervades Firedrake. Unless real function spaces play nicely with >>>>>>>>>>> communicators, Firedrake will simply crash. I'm sorry that this feature is >>>>>>>>>>> not available right now, but it actually does require significant work. >>>>>>>>>>> >>>>>>>>>>> On Wed, 10 Aug 2016 at 16:51 Anna Kalogirou < >>>>>>>>>>> A.Kalogirou(a)leeds.ac.uk> wrote: >>>>>>>>>>> >>>>>>>>>>>> Yes indeed, but I believe David would want to make it work for >>>>>>>>>>>> all cases. The merge is necessary in order to proceed with my simulations, >>>>>>>>>>>> so an update as soon as possible would be greatly appreciated. >>>>>>>>>>>> >>>>>>>>>>>> Best, Anna. >>>>>>>>>>>> >>>>>>>>>>>> >>>>>>>>>>>> On 10 Aug 2016, at 16:29, Onno Bokhove <O.Bokhove(a)leeds.ac.uk> >>>>>>>>>>>> wrote: >>>>>>>>>>>> >>>>>>>>>>>> We would already be helped for the case without MPI, I think >>>>>>>>>>>> that is correct, Anna? >>>>>>>>>>>> >>>>>>>>>>>> >>>>>>>>>>>> *From:* firedrake-bounces(a)imperial.ac.uk < >>>>>>>>>>>> firedrake-bounces(a)imperial.ac.uk> on behalf of David Ham < >>>>>>>>>>>> David.Ham(a)imperial.ac.uk> >>>>>>>>>>>> >>>>>>>>>>>> *Sent:* Wednesday, August 10, 2016 10:01 AM >>>>>>>>>>>> >>>>>>>>>>>> >>>>>>>>>>>> *To:* firedrake >>>>>>>>>>>> *Subject:* Re: [firedrake] Mixed system >>>>>>>>>>>> >>>>>>>>>>>> The short version is it depends on when I find some time to >>>>>>>>>>>> work on it. It is not simply the case that we have new features lying >>>>>>>>>>>> around waiting to be merged. The current proximate issue is to make the >>>>>>>>>>>> real spaces work properly with MPI communicators. I can't guarantee that >>>>>>>>>>>> will be the last issue. >>>>>>>>>>>> >>>>>>>>>>>> >>>>>>>>>>>> On Tue, 9 Aug 2016 at 17:37 Onno Bokhove <O.Bokhove(a)leeds.ac.uk> >>>>>>>>>>>> wrote: >>>>>>>>>>>> >>>>>>>>>>>>> It is possible to merge the branch for the mixed system >>>>>>>>>>>>> relatively soon? >>>>>>>>>>>>> Anna's results are relevant for the grant application Colin >>>>>>>>>>>>> and I aim to submit in Sept. on further wave-ship modelling. >>>>>>>>>>>>> Sorry for asking! >>>>>>>>>>>>> ------------------------------ >>>>>>>>>>>>> *From:* firedrake-bounces(a)imperial.ac.uk < >>>>>>>>>>>>> firedrake-bounces(a)imperial.ac.uk> on behalf of David Ham < >>>>>>>>>>>>> David.Ham(a)imperial.ac.uk> >>>>>>>>>>>>> *Sent:* Friday, August 5, 2016 11:37:23 AM >>>>>>>>>>>>> *To:* firedrake >>>>>>>>>>>>> *Subject:* Re: [firedrake] Mixed system >>>>>>>>>>>>> >>>>>>>>>>>>> >>>>>>>>>>>>> >>>>>>>>>>>>> On Fri, 5 Aug 2016 at 11:31 Anna Kalogirou < >>>>>>>>>>>>> a.kalogirou(a)leeds.ac.uk> wrote: >>>>>>>>>>>>> >>>>>>>>>>>>>> Hi Kramer, >>>>>>>>>>>>>> >>>>>>>>>>>>>> Thank you for the reply. >>>>>>>>>>>>>> >>>>>>>>>>>>>> The only problem I see with the solution you suggest, is that >>>>>>>>>>>>>> both I and >>>>>>>>>>>>>> lambda will be trial functions, so both terms should be on >>>>>>>>>>>>>> the bilinear >>>>>>>>>>>>>> part of the variational form. >>>>>>>>>>>>>> >>>>>>>>>>>>> >>>>>>>>>>>>> Not quite. I and lambda together form a Function in a mixed >>>>>>>>>>>>> finite element space. You solve for both of them together, just like >>>>>>>>>>>>> pressure and velocity in a fluid simulation. >>>>>>>>>>>>> >>>>>>>>>>>>> >>>>>>>>>>>>>> >>>>>>>>>>>>>> I could, of course, define a residual and solve everything >>>>>>>>>>>>>> using >>>>>>>>>>>>>> nonlinear solvers, but I don't know how efficient that would >>>>>>>>>>>>>> be. >>>>>>>>>>>>>> >>>>>>>>>>>>> >>>>>>>>>>>>> I think you are confused. None of the relevant issues change >>>>>>>>>>>>> in the nonlinear case. After all, the nonlinear solver will differentiate >>>>>>>>>>>>> the residual to create a bilinear form anyway. >>>>>>>>>>>>> >>>>>>>>>>>>> Incidentally, there is no particular reason to believe that >>>>>>>>>>>>> writing the problem in nonlinear form would be inefficient. >>>>>>>>>>>>> >>>>>>>>>>>>> Regards, >>>>>>>>>>>>> >>>>>>>>>>>>> David >>>>>>>>>>>>> >>>>>>>>>>>>> >>>>>>>>>>>>>> >>>>>>>>>>>>>> Best, >>>>>>>>>>>>>> >>>>>>>>>>>>>> Anna. >>>>>>>>>>>>>> >>>>>>>>>>>>>> Dr Anna Kalogirou >>>>>>>>>>>>>> Research Fellow >>>>>>>>>>>>>> School of Mathematics >>>>>>>>>>>>>> University of Leeds >>>>>>>>>>>>>> >>>>>>>>>>>>>> http://www1.maths.leeds.ac.uk/~matak/ >>>>>>>>>>>>>> >>>>>>>>>>>>>> On 04/08/16 17:29, Stephan Kramer wrote: >>>>>>>>>>>>>> > On 27/07/16 11:01, Anna Kalogirou wrote: >>>>>>>>>>>>>> >> Dear all, >>>>>>>>>>>>>> >> >>>>>>>>>>>>>> >> I am trying to solve the system found in the attached pdf >>>>>>>>>>>>>> as a mixed >>>>>>>>>>>>>> >> system. I want to solve the same >>>>>>>>>>>>>> >> problem as the one found here >>>>>>>>>>>>>> >> < >>>>>>>>>>>>>> https://bitbucket.org/annakalog/buoy2d/src/14334d3c20b9f10ed7c1246cde9e3cb6… >>>>>>>>>>>>>> >, >>>>>>>>>>>>>> >> >>>>>>>>>>>>>> >> but with the use of Schur complements and not linear >>>>>>>>>>>>>> algebra. >>>>>>>>>>>>>> >> >>>>>>>>>>>>>> >> In a previous discussion about this problem, Lawrence >>>>>>>>>>>>>> mentioned that the >>>>>>>>>>>>>> >> test function for integral(lambda*Theta(x-Lp)dx) needs to >>>>>>>>>>>>>> be considered >>>>>>>>>>>>>> >> as coming from the real space of constant functions. Could >>>>>>>>>>>>>> you please >>>>>>>>>>>>>> >> elaborate on this? >>>>>>>>>>>>>> > >>>>>>>>>>>>>> > My guess at what Lawrence means is (it's how I would >>>>>>>>>>>>>> probably approach >>>>>>>>>>>>>> > it), is that you add an additional equation that says >>>>>>>>>>>>>> > >>>>>>>>>>>>>> > (1) I = integral(lambda*Theta(x-Lp)dx) >>>>>>>>>>>>>> > >>>>>>>>>>>>>> > and then substitute I in your third equation. Equation (1) >>>>>>>>>>>>>> can >>>>>>>>>>>>>> > be turned into a finite element weak formulation by >>>>>>>>>>>>>> considering the >>>>>>>>>>>>>> > function space of functions that are constant over the >>>>>>>>>>>>>> entire domain >>>>>>>>>>>>>> > This function space is just 1-dimensional and therefore >>>>>>>>>>>>>> equivalent to >>>>>>>>>>>>>> > the space of reals R and hence we denote this function >>>>>>>>>>>>>> space of >>>>>>>>>>>>>> > constants simply as R. Then we can consider I to be a trial >>>>>>>>>>>>>> function >>>>>>>>>>>>>> > in R, and with a test function v4 in R, we can write: >>>>>>>>>>>>>> > >>>>>>>>>>>>>> > v4*lambda*Theta(x-Lp)*dx == v4*I*dx/area(domain) >>>>>>>>>>>>>> > >>>>>>>>>>>>>> > This is a proper finite element equation that you could >>>>>>>>>>>>>> solve in >>>>>>>>>>>>>> > conjuction with the other equations. I believe however that >>>>>>>>>>>>>> this >>>>>>>>>>>>>> > function space R is not currently implemented in Firedrake. >>>>>>>>>>>>>> It is >>>>>>>>>>>>>> > available in fenics, and is on the wishlist. >>>>>>>>>>>>>> > >>>>>>>>>>>>>> > Cheers >>>>>>>>>>>>>> > Stephan >>>>>>>>>>>>>> > >>>>>>>>>>>>>> > _______________________________________________ >>>>>>>>>>>>>> > firedrake mailing list >>>>>>>>>>>>>> > firedrake(a)imperial.ac.uk >>>>>>>>>>>>>> > https://mailman.ic.ac.uk/mailman/listinfo/firedrake >>>>>>>>>>>>>> >>>>>>>>>>>>>> >>>>>>>>>>>>>> _______________________________________________ >>>>>>>>>>>>>> firedrake mailing list >>>>>>>>>>>>>> firedrake(a)imperial.ac.uk >>>>>>>>>>>>>> https://mailman.ic.ac.uk/mailman/listinfo/firedrake >>>>>>>>>>>>>> >>>>>>>>>>>>> _______________________________________________ >>>>>>>>>>>> firedrake mailing list >>>>>>>>>>>> firedrake(a)imperial.ac.uk >>>>>>>>>>>> https://mailman.ic.ac.uk/mailman/listinfo/firedrake >>>>>>>>>>>> >>>>>>>>>>>> _______________________________________________ >>>>>>>>>> firedrake mailing list >>>>>>>>>> firedrake(a)imperial.ac.uk >>>>>>>>>> https://mailman.ic.ac.uk/mailman/listinfo/firedrake >>>>>>>>>> >>>>>>>>>> _______________________________________________ >>>>>>>>>> firedrake mailing list >>>>>>>>>> firedrake(a)imperial.ac.uk >>>>>>>>>> https://mailman.ic.ac.uk/mailman/listinfo/firedrake >>>>>>>>>> >>>>>>>>>> _______________________________________________ >>>>>>>>> firedrake mailing list >>>>>>>>> firedrake(a)imperial.ac.uk >>>>>>>>> https://mailman.ic.ac.uk/mailman/listinfo/firedrake >>>>>>>>> >>>>>>>>> _______________________________________________ >>>>>>>> firedrake mailing list >>>>>>>> firedrake(a)imperial.ac.uk >>>>>>>> https://mailman.ic.ac.uk/mailman/listinfo/firedrake >>>>>>>> >>>>>>>> _______________________________________________ >>>>> 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 >>>> >>>> _______________________________________________ >>> firedrake mailing list >>> firedrake(a)imperial.ac.uk >>> https://mailman.ic.ac.uk/mailman/listinfo/firedrake >>> >>> >>> _______________________________________________ >> firedrake mailing list >> firedrake(a)imperial.ac.uk >> https://mailman.ic.ac.uk/mailman/listinfo/firedrake >> >> >> _______________________________________________ > firedrake mailing list > firedrake(a)imperial.ac.uk > https://mailman.ic.ac.uk/mailman/listinfo/firedrake > > >

3 2

Re: [firedrake] DG for advection-diffusion
by Andrew McRae 05 Oct '16

05 Oct '16

Is yseed an integer? If so, 1/yseed presumably gives 0. On 5 October 2016 at 11:41, Justin Chang <jychang48(a)gmail.com> wrote: > Correction, h_avg = Constant(1/yseed) fails regardless if quad or > triangle. Using h_avg = 0.5*(h('+')+h('-')) seems to be the only way to > have a working code, but how do I get this equivalent for quadilateral > meshes? > > Attached are the DG problems, one with triangular elements and one with > quadrilateral elements > > On Wed, Oct 5, 2016 at 5:27 AM, Andrew McRae <A.T.T.McRae(a)bath.ac.uk> > wrote: > >> >> On 5 October 2016 at 11:25, Justin Chang <jychang48(a)gmail.com> wrote: >> >>> So I should have noted this earlier, but the code fails only with >>> quadilateral meshes. The code works fine on triangular meshes. So if I did >>> this: >>> >>> h = CellSize(mesh) >>> h_avg = 0.5*(h('+')+h('-')) >>> >>> on a quadrilateral mesh, i get this error: >>> >>> UFL:WARNING Only know how to compute the circumradius of an affine cell. >>> UFL:WARNING Only know how to compute the circumradius of an affine cell. >>> Traceback (most recent call last): >>> File "DG_test.py", line 174, in <module> >>> solve(a_D==L_D,u0_A,bcs=bcs,solver_parameters=OS_parameters, >>> options_prefix="D_") >>> File "/home/justin/Software/firedrake/src/firedrake/firedrake/solving.py", >>> line 119, in solve >>> _solve_varproblem(*args, **kwargs) >>> File "/home/justin/Software/firedrake/src/firedrake/firedrake/solving.py", >>> line 146, in _solve_varproblem >>> appctx=appctx) >>> File "/home/justin/Software/firedrake/src/firedrake/firedrake/variational_solver.py", >>> line 251, in __init__ >>> super(LinearVariationalSolver, self).__init__(*args, **kwargs) >>> File "/home/justin/Software/firedrake/src/firedrake/firedrake/variational_solver.py", >>> line 126, in __init__ >>> appctx=appctx) >>> File "/home/justin/Software/firedrake/src/firedrake/firedrake/solving_utils.py", >>> line 178, in __init__ >>> appctx=appctx) >>> File "/home/justin/Software/firedrake/src/firedrake/firedrake/assemble.py", >>> line 81, in assemble >>> inverse=inverse, mat_type=mat_type, appctx=appctx) >>> File "<decorator-gen-291>", line 2, in _assemble >>> File "/home/justin/Software/firedrake/src/firedrake/firedrake/utils.py", >>> line 62, in wrapper >>> return f(*args, **kwargs) >>> File "/home/justin/Software/firedrake/src/firedrake/firedrake/assemble.py", >>> line 120, in _assemble >>> inverse=inverse) >>> File "/home/justin/Software/firedrake/src/firedrake/firedrake/tsfc_interface.py", >>> line 189, in compile_form >>> number_map).kernels >>> File "/home/justin/Software/firedrake/src/PyOP2/pyop2/caching.py", >>> line 198, in __new__ >>> obj = make_obj() >>> File "/home/justin/Software/firedrake/src/PyOP2/pyop2/caching.py", >>> line 188, in make_obj >>> obj.__init__(*args, **kwargs) >>> File "/home/justin/Software/firedrake/src/firedrake/firedrake/tsfc_interface.py", >>> line 110, in __init__ >>> tree = tsfc_compile_form(form, prefix=name, parameters=parameters) >>> File "/home/justin/Software/firedrake/src/tsfc/tsfc/driver.py", line >>> 43, in compile_form >>> kernels.append(compile_integral(integral_data, fd, prefix, >>> parameters)) >>> File "/home/justin/Software/firedrake/src/tsfc/tsfc/driver.py", line >>> 217, in compile_integral >>> facetarea=facetarea) >>> File "/home/justin/Software/firedrake/src/tsfc/tsfc/fem.py", line >>> 405, in compile_ufl >>> return map_expr_dags(translator, expressions) >>> File "/home/justin/Software/firedrake/src/ufl/ufl/corealg/map_dag.py", >>> line 84, in map_expr_dags >>> r = handlers[v._ufl_typecode_](v) >>> File "/home/justin/Software/firedrake/src/tsfc/tsfc/ufl_utils.py", >>> line 107, in _modified_terminal >>> return self.modified_terminal(o) >>> File "/home/justin/Software/firedrake/src/tsfc/tsfc/fem.py", line >>> 254, in modified_terminal >>> return translate(mt.terminal, mt, self.parameters) >>> File "/home/justin/Software/firedrake/local/lib/python2.7/site-packages/singledispatch.py", >>> line 210, in wrapper >>> return dispatch(args[0].__class__)(*args, **kw) >>> File "/home/justin/Software/firedrake/src/tsfc/tsfc/fem.py", line >>> 312, in translate_geometricquantity >>> return geometric.translate(terminal, mt, params) >>> File "/home/justin/Software/firedrake/local/lib/python2.7/site-packages/singledispatch.py", >>> line 210, in wrapper >>> return dispatch(args[0].__class__)(*args, **kw) >>> File "/home/justin/Software/firedrake/src/tsfc/tsfc/geometric.py", >>> line 29, in translate >>> raise AssertionError("Cannot handle geometric quantity type: %s" % >>> type(terminal)) >>> AssertionError: Cannot handle geometric quantity type: <class >>> 'ufl.geometry.Circumradius'> >>> >>> So I resorted to: >>> >>> h_avg = Constant(1/yseed) >>> >>> which also works for triangular meshes but not work quadrilateral meshes. >>> >> >> What do you mean by 'not work'? >> >> >>> >>> Is there a more appropriate way to obtain the cell size for quads? >>> >>> Justin >>> >>> On Wed, Oct 5, 2016 at 4:41 AM, David Ham <David.Ham(a)imperial.ac.uk> >>> wrote: >>> >>>> Hi Justin, >>>> >>>> A couple of things to try: >>>> >>>> 1. Work with the weak BCs. The strong BC formulation is rather dubious >>>> from a maths perspective (not that we expect it to be wrong, but still...) >>>> >>>> 2. Cut down to single material in the first instance. >>>> >>>> 3. Drop back to 1D so you can look at actual matrices, and convince >>>> yourself that the matrix you are producing has the right invariants. EG >>>> positive definite, known couplings correct etc. >>>> >>>> On Wed, 5 Oct 2016 at 10:14 Justin Chang <jychang48(a)gmail.com> wrote: >>>> >>>>> Sorry, here's a simpler version. Still has the same errors regardless >>>>> of whether bcs are strongly or weakly imposed. Direct solver (-D_ksp_type >>>>> preonly -D_pc_type lu) gives a complete nonsensical answer. Run problem as: >>>>> >>>>> python DG_test.py 100 50 >>>>> >>>>> Any help appreciated, thanks! >>>>> >>>>> Justin >>>>> >>>>> On Wed, Oct 5, 2016 at 3:28 AM, Andrew McRae <A.T.T.McRae(a)bath.ac.uk> >>>>> wrote: >>>>> >>>>> Hi Justin, >>>>> >>>>> I don't think it's a good use of time for us to debug your full >>>>> program (which runs an external python file at line 74, for example!) >>>>> >>>>> I suggest you strip down the program to its absolute core. If even >>>>> this doesn't run then get back in touch with a minimal failing example. >>>>> >>>>> Thanks, >>>>> Andrew >>>>> >>>>> On 5 October 2016 at 09:22, Justin Chang <jychang48(a)gmail.com> wrote: >>>>> >>>>> Thanks, so the problem was that the diffusivity tensor was not inside >>>>> the avg() terms. That problem is fixed. >>>>> >>>>> However, when the code was modified as such, i get this error: >>>>> >>>>> Traceback (most recent call last): >>>>> File "2D_OS_boundary_ex1.py", line 313, in <module> >>>>> u0 = ADsolve(solver_D,tao_D,opt_D,u0) >>>>> File "2D_OS_boundary_ex1.py", line 292, in ADsolve >>>>> solver_O.solve(u0_O,b) >>>>> File "/home/justin/Software/firedrake/src/firedrake/firedrake/linear_solver.py", >>>>> line 141, in solve >>>>> raise ConvergenceError("LinearSolver failed to converge after %d >>>>> iterations with reason: %s", self.ksp.getIterationNumber(), >>>>> solving_utils.KSPReasons[r]) >>>>> firedrake.exceptions.ConvergenceError: ('LinearSolver failed to >>>>> converge after %d iterations with reason: %s', 0, 'DIVERGED_NANORINF') >>>>> >>>>> Is something wrong with my weak formulation? I even tried switching >>>>> back to weakly applied Boundaries but I still get the same error. >>>>> >>>>> Thanks, >>>>> Justin >>>>> >>>>> On Wed, Oct 5, 2016 at 3:05 AM, Andrew McRae <A.T.T.McRae(a)bath.ac.uk> >>>>> wrote: >>>>> >>>>> Because all TrialFunctions, TestFunctions and Functions in interior >>>>> facet integrals need to include restrictions (or use operators such as jump >>>>> and avg) >>>>> >>>>> On 5 October 2016 at 08:32, Justin Chang <jychang48(a)gmail.com> wrote: >>>>> >>>>> Okay thanks, >>>>> >>>>> I am attempting to do the diffusion part only at the moment under DG, >>>>> but I am getting these errors: >>>>> >>>>> $ python 2D_OS_boundary_ex1.py 100 50 0 0 0 -D_pc_type bjacobi >>>>> UFL:WARNING Discontinuous Lagrange element requested on quadrilateral, >>>>> creating DQ element. >>>>> UFL:WARNING Discontinuous Lagrange element requested on quadrilateral, >>>>> creating DQ element. >>>>> UFL:WARNING Discontinuous Lagrange element requested on quadrilateral, >>>>> creating DQ element. >>>>> UFL:ERROR Discontinuous type Coefficient must be restricted. >>>>> Traceback (most recent call last): >>>>> File "2D_OS_boundary_ex1.py", line 255, in <module> >>>>> A_D = assemble(a_D,bcs=bc_D) >>>>> File "/home/justin/Software/firedrake/src/firedrake/firedrake/assemble.py", >>>>> line 81, in assemble >>>>> inverse=inverse, mat_type=mat_type, appctx=appctx) >>>>> File "<decorator-gen-291>", line 2, in _assemble >>>>> File "/home/justin/Software/firedrake/src/firedrake/firedrake/utils.py", >>>>> line 62, in wrapper >>>>> return f(*args, **kwargs) >>>>> File "/home/justin/Software/firedrake/src/firedrake/firedrake/assemble.py", >>>>> line 120, in _assemble >>>>> inverse=inverse) >>>>> File "/home/justin/Software/firedrake/src/firedrake/firedrake/tsfc_interface.py", >>>>> line 189, in compile_form >>>>> number_map).kernels >>>>> File "/home/justin/Software/firedrake/src/PyOP2/pyop2/caching.py", >>>>> line 198, in __new__ >>>>> obj = make_obj() >>>>> File "/home/justin/Software/firedrake/src/PyOP2/pyop2/caching.py", >>>>> line 188, in make_obj >>>>> obj.__init__(*args, **kwargs) >>>>> File "/home/justin/Software/firedrake/src/firedrake/firedrake/tsfc_interface.py", >>>>> line 110, in __init__ >>>>> tree = tsfc_compile_form(form, prefix=name, parameters=parameters) >>>>> File "/home/justin/Software/firedrake/src/tsfc/tsfc/driver.py", >>>>> line 36, in compile_form >>>>> fd = ufl_utils.compute_form_data(form) >>>>> File "/home/justin/Software/firedrake/src/tsfc/tsfc/ufl_utils.py", >>>>> line 47, in compute_form_data >>>>> do_estimate_degrees=do_estimate_degrees, >>>>> File "/home/justin/Software/firedrake/src/ufl/ufl/algorithms/compute_form_data.py", >>>>> line 305, in compute_form_data >>>>> form = apply_restrictions(form) >>>>> File "/home/justin/Software/firedrake/src/ufl/ufl/algorithms/apply_restrictions.py", >>>>> line 211, in apply_restrictions >>>>> only_integral_type=integral_types) >>>>> File "/home/justin/Software/firedrake/src/ufl/ufl/algorithms/map_integrands.py", >>>>> line 58, in map_integrand_dags >>>>> form, only_integral_type) >>>>> File "/home/justin/Software/firedrake/src/ufl/ufl/algorithms/map_integrands.py", >>>>> line 39, in map_integrands >>>>> for itg in form.integrals()] >>>>> File "/home/justin/Software/firedrake/src/ufl/ufl/algorithms/map_integrands.py", >>>>> line 46, in map_integrands >>>>> return itg.reconstruct(function(itg.integrand())) >>>>> File "/home/justin/Software/firedrake/src/ufl/ufl/algorithms/map_integrands.py", >>>>> line 57, in <lambda> >>>>> return map_integrands(lambda expr: map_expr_dag(function, expr, >>>>> compress), >>>>> File "/home/justin/Software/firedrake/src/ufl/ufl/corealg/map_dag.py", >>>>> line 37, in map_expr_dag >>>>> result, = map_expr_dags(function, [expression], compress=compress) >>>>> File "/home/justin/Software/firedrake/src/ufl/ufl/corealg/map_dag.py", >>>>> line 84, in map_expr_dags >>>>> r = handlers[v._ufl_typecode_](v) >>>>> File "/home/justin/Software/firedrake/src/ufl/ufl/algorithms/apply_restrictions.py", >>>>> line 102, in reference_value >>>>> g = self(f) >>>>> File "/home/justin/Software/firedrake/src/ufl/ufl/corealg/multifunction.py", >>>>> line 95, in __call__ >>>>> return self._handlers[o._ufl_typecode_](o, *args) >>>>> File "/home/justin/Software/firedrake/src/ufl/ufl/algorithms/apply_restrictions.py", >>>>> line 129, in coefficient >>>>> return self._require_restriction(o) >>>>> File "/home/justin/Software/firedrake/src/ufl/ufl/algorithms/apply_restrictions.py", >>>>> line 59, in _require_restriction >>>>> "Discontinuous type %s must be restricted." % >>>>> o._ufl_class_.__name__) >>>>> File "/home/justin/Software/firedrake/src/ufl/ufl/assertions.py", >>>>> line 51, in ufl_assert >>>>> error(*message) >>>>> File "/home/justin/Software/firedrake/src/ufl/ufl/log.py", line >>>>> 167, in error >>>>> raise self._exception_type(self._format_raw(*message)) >>>>> ufl.log.UFLException: Discontinuous type Coefficient must be >>>>> restricted. >>>>> >>>>> I can't seem to figure out why. Attached is the code. Run as: >>>>> >>>>> python 2D_OS_boundary_ex1.py 100 50 0 0 0 >>>>> >>>>> Thanks, >>>>> Justin >>>>> >>>>> On Thu, Sep 29, 2016 at 4:49 AM, David Ham <David.Ham(a)imperial.ac.uk> >>>>> wrote: >>>>> >>>>> >>>>> >>>>> On Thu, 29 Sep 2016 at 10:34 Justin Chang <jychang48(a)gmail.com> wrote: >>>>> >>>>> Hi all, >>>>> >>>>> I have written here the Discontinuous Galerkin method for 1D >>>>> advection-diffusion equations using backward euler time discretization: >>>>> >>>>> mesh = UnitIntervalMesh(100) >>>>> V = FunctionSpace(mesh, 'DG', 1) >>>>> u = TrialFunction(V) >>>>> v = TestFunction(V) >>>>> >>>>> u0 = Function(V, name="initial") >>>>> sol = Function(V, name="solution") >>>>> >>>>> # Diffusivity >>>>> D = Constant(0.005) >>>>> >>>>> # Velocity >>>>> vel = as_vector((Constant(1.0),)) >>>>> >>>>> # Forcing function >>>>> f = interpolate(Expression(....),V) >>>>> >>>>> # Penalty >>>>> k = 1 >>>>> d = 1 >>>>> gamma = Constant((k+1)*(k+d)/d)*FacetArea(mesh)/avg(CellVolume(mesh)) >>>>> # based on the formula in http://webpages.sdsmt.edu/~ksh >>>>> ahbaz/shahbazi_05.pdf >>>>> >>>>> # Other parameters >>>>> n = FacetNormal(mesh) >>>>> h = CellSize(mesh) >>>>> vn = 0.5*(dot(vel,n) + abs(dot(vel,n))) >>>>> >>>>> # Weak form >>>>> a = v*u/dt*dx - dot(grad(v),vel*u)*dx + dot(jump(v),vn('+')*u('+')-vn('-')*u('-'))*dS >>>>> + dot(v, vn*u)*ds \ >>>>> + dot(grad(v),D*grad(u))*dx - dot(jump(v,n),avg(D*grad(u)))*dS - >>>>> dot(avg(D*grad(v)),jump(u,n))*dS \ >>>>> + gamma*dot(jump(v,n),jump(u,n))*dS >>>>> L = v*f*dx + v*u0/dt*dx >>>>> >>>>> >>>>> ... >>>>> >>>>> I have a couple questions: >>>>> >>>>> 1) If say I am solving for the transport of 4 chemical species, >>>>> Normally I would create a mixed function space like V_mixed = V*V*V*V, >>>>> split u into u1, u2, u3, and u4, write bilinear and linear forms for each >>>>> of these individual components and add them up. That is, a = a1 + a2 + a3 + >>>>> a4, L = L1 + L2 + L3 + L4. This can become quite a lot of code if I end up >>>>> having to solve ~20 chemical species eventually. The catch is that the >>>>> velocity and diffusivity will be the same for all species but the forcing >>>>> function and boundary conditions are going to be different. Is there a more >>>>> efficient of writing the bilinear and linear form instead of having to copy >>>>> and paste the a = ... and L = ... multiple times and changing the >>>>> trial/test functions? >>>>> >>>>> >>>>> a) If you are independently advecting lots of species (ie no reaction) >>>>> then you probably want to just solve the equations separately, rather than >>>>> making a massive mixed function space. >>>>> >>>>> b) You can avoid lots of retyping by just using normal Python >>>>> functions to construct forms. You could, for example, write a function >>>>> which takes a test and trial function, and maybe some BC information, and >>>>> returns the appropriate forms. >>>>> >>>>> >>>>> >>>>> 2) I see that in FEniCS's undocumented dg-advection-diffusion example, >>>>> they enforce the DirichletBC's strongly. Is it possible to enforce them >>>>> weakly like the SIPG or DG upwinding techniques for the pure diffusion or >>>>> advection equations? >>>>> >>>>> >>>>> Yes, just substitute boundary conditions for surface integrals in the >>>>> usual way. >>>>> >>>>> >>>>> >>>>> >>>>> Thanks, >>>>> Justin >>>>> >>>>> >>>>> >>>>> _______________________________________________ >>>>> firedrake mailing list >>>>> firedrake(a)imperial.ac.uk >>>>> https://mailman.ic.ac.uk/mailman/listinfo/firedrake >>>>> >>>>> >>>>> _______________________________________________ >>>>> firedrake mailing list >>>>> firedrake(a)imperial.ac.uk >>>>> https://mailman.ic.ac.uk/mailman/listinfo/firedrake >>>>> >>>>> >>>>> _______________________________________________ >>>>> firedrake mailing list >>>>> firedrake(a)imperial.ac.uk >>>>> https://mailman.ic.ac.uk/mailman/listinfo/firedrake >>>>> >>>>> >>>> _______________________________________________ >>>> firedrake mailing list >>>> firedrake(a)imperial.ac.uk >>>> https://mailman.ic.ac.uk/mailman/listinfo/firedrake >>>> >>>> >>> >> >> _______________________________________________ >> firedrake mailing list >> firedrake(a)imperial.ac.uk >> https://mailman.ic.ac.uk/mailman/listinfo/firedrake >> >> >

2 2

Re: [firedrake] DG for advection-diffusion
by Andrew McRae 05 Oct '16

05 Oct '16

On 5 October 2016 at 11:25, Justin Chang <jychang48(a)gmail.com> wrote: > So I should have noted this earlier, but the code fails only with > quadilateral meshes. The code works fine on triangular meshes. So if I did > this: > > h = CellSize(mesh) > h_avg = 0.5*(h('+')+h('-')) > > on a quadrilateral mesh, i get this error: > > UFL:WARNING Only know how to compute the circumradius of an affine cell. > UFL:WARNING Only know how to compute the circumradius of an affine cell. > Traceback (most recent call last): > File "DG_test.py", line 174, in <module> > solve(a_D==L_D,u0_A,bcs=bcs,solver_parameters=OS_ > parameters,options_prefix="D_") > File "/home/justin/Software/firedrake/src/firedrake/firedrake/solving.py", > line 119, in solve > _solve_varproblem(*args, **kwargs) > File "/home/justin/Software/firedrake/src/firedrake/firedrake/solving.py", > line 146, in _solve_varproblem > appctx=appctx) > File "/home/justin/Software/firedrake/src/firedrake/ > firedrake/variational_solver.py", line 251, in __init__ > super(LinearVariationalSolver, self).__init__(*args, **kwargs) > File "/home/justin/Software/firedrake/src/firedrake/ > firedrake/variational_solver.py", line 126, in __init__ > appctx=appctx) > File "/home/justin/Software/firedrake/src/firedrake/firedrake/solving_utils.py", > line 178, in __init__ > appctx=appctx) > File "/home/justin/Software/firedrake/src/firedrake/firedrake/assemble.py", > line 81, in assemble > inverse=inverse, mat_type=mat_type, appctx=appctx) > File "<decorator-gen-291>", line 2, in _assemble > File "/home/justin/Software/firedrake/src/firedrake/firedrake/utils.py", > line 62, in wrapper > return f(*args, **kwargs) > File "/home/justin/Software/firedrake/src/firedrake/firedrake/assemble.py", > line 120, in _assemble > inverse=inverse) > File "/home/justin/Software/firedrake/src/firedrake/firedrake/tsfc_interface.py", > line 189, in compile_form > number_map).kernels > File "/home/justin/Software/firedrake/src/PyOP2/pyop2/caching.py", line > 198, in __new__ > obj = make_obj() > File "/home/justin/Software/firedrake/src/PyOP2/pyop2/caching.py", line > 188, in make_obj > obj.__init__(*args, **kwargs) > File "/home/justin/Software/firedrake/src/firedrake/firedrake/tsfc_interface.py", > line 110, in __init__ > tree = tsfc_compile_form(form, prefix=name, parameters=parameters) > File "/home/justin/Software/firedrake/src/tsfc/tsfc/driver.py", line > 43, in compile_form > kernels.append(compile_integral(integral_data, fd, prefix, > parameters)) > File "/home/justin/Software/firedrake/src/tsfc/tsfc/driver.py", line > 217, in compile_integral > facetarea=facetarea) > File "/home/justin/Software/firedrake/src/tsfc/tsfc/fem.py", line 405, > in compile_ufl > return map_expr_dags(translator, expressions) > File "/home/justin/Software/firedrake/src/ufl/ufl/corealg/map_dag.py", > line 84, in map_expr_dags > r = handlers[v._ufl_typecode_](v) > File "/home/justin/Software/firedrake/src/tsfc/tsfc/ufl_utils.py", line > 107, in _modified_terminal > return self.modified_terminal(o) > File "/home/justin/Software/firedrake/src/tsfc/tsfc/fem.py", line 254, > in modified_terminal > return translate(mt.terminal, mt, self.parameters) > File "/home/justin/Software/firedrake/local/lib/python2.7/ > site-packages/singledispatch.py", line 210, in wrapper > return dispatch(args[0].__class__)(*args, **kw) > File "/home/justin/Software/firedrake/src/tsfc/tsfc/fem.py", line 312, > in translate_geometricquantity > return geometric.translate(terminal, mt, params) > File "/home/justin/Software/firedrake/local/lib/python2.7/ > site-packages/singledispatch.py", line 210, in wrapper > return dispatch(args[0].__class__)(*args, **kw) > File "/home/justin/Software/firedrake/src/tsfc/tsfc/geometric.py", line > 29, in translate > raise AssertionError("Cannot handle geometric quantity type: %s" % > type(terminal)) > AssertionError: Cannot handle geometric quantity type: <class > 'ufl.geometry.Circumradius'> > > So I resorted to: > > h_avg = Constant(1/yseed) > > which also works for triangular meshes but not work quadrilateral meshes. > What do you mean by 'not work'? > > Is there a more appropriate way to obtain the cell size for quads? > > Justin > > On Wed, Oct 5, 2016 at 4:41 AM, David Ham <David.Ham(a)imperial.ac.uk> > wrote: > >> Hi Justin, >> >> A couple of things to try: >> >> 1. Work with the weak BCs. The strong BC formulation is rather dubious >> from a maths perspective (not that we expect it to be wrong, but still...) >> >> 2. Cut down to single material in the first instance. >> >> 3. Drop back to 1D so you can look at actual matrices, and convince >> yourself that the matrix you are producing has the right invariants. EG >> positive definite, known couplings correct etc. >> >> On Wed, 5 Oct 2016 at 10:14 Justin Chang <jychang48(a)gmail.com> wrote: >> >>> Sorry, here's a simpler version. Still has the same errors regardless of >>> whether bcs are strongly or weakly imposed. Direct solver (-D_ksp_type >>> preonly -D_pc_type lu) gives a complete nonsensical answer. Run problem as: >>> >>> python DG_test.py 100 50 >>> >>> Any help appreciated, thanks! >>> >>> Justin >>> >>> On Wed, Oct 5, 2016 at 3:28 AM, Andrew McRae <A.T.T.McRae(a)bath.ac.uk> >>> wrote: >>> >>> Hi Justin, >>> >>> I don't think it's a good use of time for us to debug your full program >>> (which runs an external python file at line 74, for example!) >>> >>> I suggest you strip down the program to its absolute core. If even this >>> doesn't run then get back in touch with a minimal failing example. >>> >>> Thanks, >>> Andrew >>> >>> On 5 October 2016 at 09:22, Justin Chang <jychang48(a)gmail.com> wrote: >>> >>> Thanks, so the problem was that the diffusivity tensor was not inside >>> the avg() terms. That problem is fixed. >>> >>> However, when the code was modified as such, i get this error: >>> >>> Traceback (most recent call last): >>> File "2D_OS_boundary_ex1.py", line 313, in <module> >>> u0 = ADsolve(solver_D,tao_D,opt_D,u0) >>> File "2D_OS_boundary_ex1.py", line 292, in ADsolve >>> solver_O.solve(u0_O,b) >>> File "/home/justin/Software/firedrake/src/firedrake/firedrake/linear_solver.py", >>> line 141, in solve >>> raise ConvergenceError("LinearSolver failed to converge after %d >>> iterations with reason: %s", self.ksp.getIterationNumber(), >>> solving_utils.KSPReasons[r]) >>> firedrake.exceptions.ConvergenceError: ('LinearSolver failed to >>> converge after %d iterations with reason: %s', 0, 'DIVERGED_NANORINF') >>> >>> Is something wrong with my weak formulation? I even tried switching back >>> to weakly applied Boundaries but I still get the same error. >>> >>> Thanks, >>> Justin >>> >>> On Wed, Oct 5, 2016 at 3:05 AM, Andrew McRae <A.T.T.McRae(a)bath.ac.uk> >>> wrote: >>> >>> Because all TrialFunctions, TestFunctions and Functions in interior >>> facet integrals need to include restrictions (or use operators such as jump >>> and avg) >>> >>> On 5 October 2016 at 08:32, Justin Chang <jychang48(a)gmail.com> wrote: >>> >>> Okay thanks, >>> >>> I am attempting to do the diffusion part only at the moment under DG, >>> but I am getting these errors: >>> >>> $ python 2D_OS_boundary_ex1.py 100 50 0 0 0 -D_pc_type bjacobi >>> UFL:WARNING Discontinuous Lagrange element requested on quadrilateral, >>> creating DQ element. >>> UFL:WARNING Discontinuous Lagrange element requested on quadrilateral, >>> creating DQ element. >>> UFL:WARNING Discontinuous Lagrange element requested on quadrilateral, >>> creating DQ element. >>> UFL:ERROR Discontinuous type Coefficient must be restricted. >>> Traceback (most recent call last): >>> File "2D_OS_boundary_ex1.py", line 255, in <module> >>> A_D = assemble(a_D,bcs=bc_D) >>> File "/home/justin/Software/firedrake/src/firedrake/firedrake/assemble.py", >>> line 81, in assemble >>> inverse=inverse, mat_type=mat_type, appctx=appctx) >>> File "<decorator-gen-291>", line 2, in _assemble >>> File "/home/justin/Software/firedrake/src/firedrake/firedrake/utils.py", >>> line 62, in wrapper >>> return f(*args, **kwargs) >>> File "/home/justin/Software/firedrake/src/firedrake/firedrake/assemble.py", >>> line 120, in _assemble >>> inverse=inverse) >>> File "/home/justin/Software/firedrake/src/firedrake/firedrake/tsfc_interface.py", >>> line 189, in compile_form >>> number_map).kernels >>> File "/home/justin/Software/firedrake/src/PyOP2/pyop2/caching.py", >>> line 198, in __new__ >>> obj = make_obj() >>> File "/home/justin/Software/firedrake/src/PyOP2/pyop2/caching.py", >>> line 188, in make_obj >>> obj.__init__(*args, **kwargs) >>> File "/home/justin/Software/firedrake/src/firedrake/firedrake/tsfc_interface.py", >>> line 110, in __init__ >>> tree = tsfc_compile_form(form, prefix=name, parameters=parameters) >>> File "/home/justin/Software/firedrake/src/tsfc/tsfc/driver.py", line >>> 36, in compile_form >>> fd = ufl_utils.compute_form_data(form) >>> File "/home/justin/Software/firedrake/src/tsfc/tsfc/ufl_utils.py", >>> line 47, in compute_form_data >>> do_estimate_degrees=do_estimate_degrees, >>> File "/home/justin/Software/firedrake/src/ufl/ufl/algorithms/compute_form_data.py", >>> line 305, in compute_form_data >>> form = apply_restrictions(form) >>> File "/home/justin/Software/firedrake/src/ufl/ufl/algorithms/apply_restrictions.py", >>> line 211, in apply_restrictions >>> only_integral_type=integral_types) >>> File "/home/justin/Software/firedrake/src/ufl/ufl/algorithms/map_integrands.py", >>> line 58, in map_integrand_dags >>> form, only_integral_type) >>> File "/home/justin/Software/firedrake/src/ufl/ufl/algorithms/map_integrands.py", >>> line 39, in map_integrands >>> for itg in form.integrals()] >>> File "/home/justin/Software/firedrake/src/ufl/ufl/algorithms/map_integrands.py", >>> line 46, in map_integrands >>> return itg.reconstruct(function(itg.integrand())) >>> File "/home/justin/Software/firedrake/src/ufl/ufl/algorithms/map_integrands.py", >>> line 57, in <lambda> >>> return map_integrands(lambda expr: map_expr_dag(function, expr, >>> compress), >>> File "/home/justin/Software/firedrake/src/ufl/ufl/corealg/map_dag.py", >>> line 37, in map_expr_dag >>> result, = map_expr_dags(function, [expression], compress=compress) >>> File "/home/justin/Software/firedrake/src/ufl/ufl/corealg/map_dag.py", >>> line 84, in map_expr_dags >>> r = handlers[v._ufl_typecode_](v) >>> File "/home/justin/Software/firedrake/src/ufl/ufl/algorithms/apply_restrictions.py", >>> line 102, in reference_value >>> g = self(f) >>> File "/home/justin/Software/firedrake/src/ufl/ufl/corealg/multifunction.py", >>> line 95, in __call__ >>> return self._handlers[o._ufl_typecode_](o, *args) >>> File "/home/justin/Software/firedrake/src/ufl/ufl/algorithms/apply_restrictions.py", >>> line 129, in coefficient >>> return self._require_restriction(o) >>> File "/home/justin/Software/firedrake/src/ufl/ufl/algorithms/apply_restrictions.py", >>> line 59, in _require_restriction >>> "Discontinuous type %s must be restricted." % o._ufl_class_.__name__) >>> File "/home/justin/Software/firedrake/src/ufl/ufl/assertions.py", >>> line 51, in ufl_assert >>> error(*message) >>> File "/home/justin/Software/firedrake/src/ufl/ufl/log.py", line 167, >>> in error >>> raise self._exception_type(self._format_raw(*message)) >>> ufl.log.UFLException: Discontinuous type Coefficient must be restricted. >>> >>> I can't seem to figure out why. Attached is the code. Run as: >>> >>> python 2D_OS_boundary_ex1.py 100 50 0 0 0 >>> >>> Thanks, >>> Justin >>> >>> On Thu, Sep 29, 2016 at 4:49 AM, David Ham <David.Ham(a)imperial.ac.uk> >>> wrote: >>> >>> >>> >>> On Thu, 29 Sep 2016 at 10:34 Justin Chang <jychang48(a)gmail.com> wrote: >>> >>> Hi all, >>> >>> I have written here the Discontinuous Galerkin method for 1D >>> advection-diffusion equations using backward euler time discretization: >>> >>> mesh = UnitIntervalMesh(100) >>> V = FunctionSpace(mesh, 'DG', 1) >>> u = TrialFunction(V) >>> v = TestFunction(V) >>> >>> u0 = Function(V, name="initial") >>> sol = Function(V, name="solution") >>> >>> # Diffusivity >>> D = Constant(0.005) >>> >>> # Velocity >>> vel = as_vector((Constant(1.0),)) >>> >>> # Forcing function >>> f = interpolate(Expression(....),V) >>> >>> # Penalty >>> k = 1 >>> d = 1 >>> gamma = Constant((k+1)*(k+d)/d)*FacetArea(mesh)/avg(CellVolume(mesh)) # >>> based on the formula in http://webpages.sdsmt.edu/~ksh >>> ahbaz/shahbazi_05.pdf >>> >>> # Other parameters >>> n = FacetNormal(mesh) >>> h = CellSize(mesh) >>> vn = 0.5*(dot(vel,n) + abs(dot(vel,n))) >>> >>> # Weak form >>> a = v*u/dt*dx - dot(grad(v),vel*u)*dx + dot(jump(v),vn('+')*u('+')-vn('-')*u('-'))*dS >>> + dot(v, vn*u)*ds \ >>> + dot(grad(v),D*grad(u))*dx - dot(jump(v,n),avg(D*grad(u)))*dS - >>> dot(avg(D*grad(v)),jump(u,n))*dS \ >>> + gamma*dot(jump(v,n),jump(u,n))*dS >>> L = v*f*dx + v*u0/dt*dx >>> >>> >>> ... >>> >>> I have a couple questions: >>> >>> 1) If say I am solving for the transport of 4 chemical species, Normally >>> I would create a mixed function space like V_mixed = V*V*V*V, split u into >>> u1, u2, u3, and u4, write bilinear and linear forms for each of these >>> individual components and add them up. That is, a = a1 + a2 + a3 + a4, L = >>> L1 + L2 + L3 + L4. This can become quite a lot of code if I end up having >>> to solve ~20 chemical species eventually. The catch is that the velocity >>> and diffusivity will be the same for all species but the forcing function >>> and boundary conditions are going to be different. Is there a more >>> efficient of writing the bilinear and linear form instead of having to copy >>> and paste the a = ... and L = ... multiple times and changing the >>> trial/test functions? >>> >>> >>> a) If you are independently advecting lots of species (ie no reaction) >>> then you probably want to just solve the equations separately, rather than >>> making a massive mixed function space. >>> >>> b) You can avoid lots of retyping by just using normal Python functions >>> to construct forms. You could, for example, write a function which takes a >>> test and trial function, and maybe some BC information, and returns the >>> appropriate forms. >>> >>> >>> >>> 2) I see that in FEniCS's undocumented dg-advection-diffusion example, >>> they enforce the DirichletBC's strongly. Is it possible to enforce them >>> weakly like the SIPG or DG upwinding techniques for the pure diffusion or >>> advection equations? >>> >>> >>> Yes, just substitute boundary conditions for surface integrals in the >>> usual way. >>> >>> >>> >>> >>> Thanks, >>> Justin >>> >>> >>> >>> _______________________________________________ >>> firedrake mailing list >>> firedrake(a)imperial.ac.uk >>> https://mailman.ic.ac.uk/mailman/listinfo/firedrake >>> >>> >>> _______________________________________________ >>> firedrake mailing list >>> firedrake(a)imperial.ac.uk >>> https://mailman.ic.ac.uk/mailman/listinfo/firedrake >>> >>> >>> _______________________________________________ >>> firedrake mailing list >>> firedrake(a)imperial.ac.uk >>> https://mailman.ic.ac.uk/mailman/listinfo/firedrake >>> >>> >> _______________________________________________ >> firedrake mailing list >> firedrake(a)imperial.ac.uk >> https://mailman.ic.ac.uk/mailman/listinfo/firedrake >> >> >

2 1

Re: [firedrake] DG for advection-diffusion
by David Ham 05 Oct '16

05 Oct '16

Hi Justin, A couple of things to try: 1. Work with the weak BCs. The strong BC formulation is rather dubious from a maths perspective (not that we expect it to be wrong, but still...) 2. Cut down to single material in the first instance. 3. Drop back to 1D so you can look at actual matrices, and convince yourself that the matrix you are producing has the right invariants. EG positive definite, known couplings correct etc. On Wed, 5 Oct 2016 at 10:14 Justin Chang <jychang48(a)gmail.com> wrote: > Sorry, here's a simpler version. Still has the same errors regardless of > whether bcs are strongly or weakly imposed. Direct solver (-D_ksp_type > preonly -D_pc_type lu) gives a complete nonsensical answer. Run problem as: > > python DG_test.py 100 50 > > Any help appreciated, thanks! > > Justin > > On Wed, Oct 5, 2016 at 3:28 AM, Andrew McRae <A.T.T.McRae(a)bath.ac.uk> > wrote: > > Hi Justin, > > I don't think it's a good use of time for us to debug your full program > (which runs an external python file at line 74, for example!) > > I suggest you strip down the program to its absolute core. If even this > doesn't run then get back in touch with a minimal failing example. > > Thanks, > Andrew > > On 5 October 2016 at 09:22, Justin Chang <jychang48(a)gmail.com> wrote: > > Thanks, so the problem was that the diffusivity tensor was not inside the > avg() terms. That problem is fixed. > > However, when the code was modified as such, i get this error: > > Traceback (most recent call last): > File "2D_OS_boundary_ex1.py", line 313, in <module> > u0 = ADsolve(solver_D,tao_D,opt_D,u0) > File "2D_OS_boundary_ex1.py", line 292, in ADsolve > solver_O.solve(u0_O,b) > File > "/home/justin/Software/firedrake/src/firedrake/firedrake/linear_solver.py", > line 141, in solve > raise ConvergenceError("LinearSolver failed to converge after %d > iterations with reason: %s", self.ksp.getIterationNumber(), > solving_utils.KSPReasons[r]) > firedrake.exceptions.ConvergenceError: ('LinearSolver failed to converge > after %d iterations with reason: %s', 0, 'DIVERGED_NANORINF') > > Is something wrong with my weak formulation? I even tried switching back > to weakly applied Boundaries but I still get the same error. > > Thanks, > Justin > > On Wed, Oct 5, 2016 at 3:05 AM, Andrew McRae <A.T.T.McRae(a)bath.ac.uk> > wrote: > > Because all TrialFunctions, TestFunctions and Functions in interior facet > integrals need to include restrictions (or use operators such as jump and > avg) > > On 5 October 2016 at 08:32, Justin Chang <jychang48(a)gmail.com> wrote: > > Okay thanks, > > I am attempting to do the diffusion part only at the moment under DG, but > I am getting these errors: > > $ python 2D_OS_boundary_ex1.py 100 50 0 0 0 -D_pc_type bjacobi > UFL:WARNING Discontinuous Lagrange element requested on quadrilateral, > creating DQ element. > UFL:WARNING Discontinuous Lagrange element requested on quadrilateral, > creating DQ element. > UFL:WARNING Discontinuous Lagrange element requested on quadrilateral, > creating DQ element. > UFL:ERROR Discontinuous type Coefficient must be restricted. > Traceback (most recent call last): > File "2D_OS_boundary_ex1.py", line 255, in <module> > A_D = assemble(a_D,bcs=bc_D) > File > "/home/justin/Software/firedrake/src/firedrake/firedrake/assemble.py", line > 81, in assemble > inverse=inverse, mat_type=mat_type, appctx=appctx) > File "<decorator-gen-291>", line 2, in _assemble > File "/home/justin/Software/firedrake/src/firedrake/firedrake/utils.py", > line 62, in wrapper > return f(*args, **kwargs) > File > "/home/justin/Software/firedrake/src/firedrake/firedrake/assemble.py", line > 120, in _assemble > inverse=inverse) > File > "/home/justin/Software/firedrake/src/firedrake/firedrake/tsfc_interface.py", > line 189, in compile_form > number_map).kernels > File "/home/justin/Software/firedrake/src/PyOP2/pyop2/caching.py", line > 198, in __new__ > obj = make_obj() > File "/home/justin/Software/firedrake/src/PyOP2/pyop2/caching.py", line > 188, in make_obj > obj.__init__(*args, **kwargs) > File > "/home/justin/Software/firedrake/src/firedrake/firedrake/tsfc_interface.py", > line 110, in __init__ > tree = tsfc_compile_form(form, prefix=name, parameters=parameters) > File "/home/justin/Software/firedrake/src/tsfc/tsfc/driver.py", line 36, > in compile_form > fd = ufl_utils.compute_form_data(form) > File "/home/justin/Software/firedrake/src/tsfc/tsfc/ufl_utils.py", line > 47, in compute_form_data > do_estimate_degrees=do_estimate_degrees, > File > "/home/justin/Software/firedrake/src/ufl/ufl/algorithms/compute_form_data.py", > line 305, in compute_form_data > form = apply_restrictions(form) > File > "/home/justin/Software/firedrake/src/ufl/ufl/algorithms/apply_restrictions.py", > line 211, in apply_restrictions > only_integral_type=integral_types) > File > "/home/justin/Software/firedrake/src/ufl/ufl/algorithms/map_integrands.py", > line 58, in map_integrand_dags > form, only_integral_type) > File > "/home/justin/Software/firedrake/src/ufl/ufl/algorithms/map_integrands.py", > line 39, in map_integrands > for itg in form.integrals()] > File > "/home/justin/Software/firedrake/src/ufl/ufl/algorithms/map_integrands.py", > line 46, in map_integrands > return itg.reconstruct(function(itg.integrand())) > File > "/home/justin/Software/firedrake/src/ufl/ufl/algorithms/map_integrands.py", > line 57, in <lambda> > return map_integrands(lambda expr: map_expr_dag(function, expr, > compress), > File "/home/justin/Software/firedrake/src/ufl/ufl/corealg/map_dag.py", > line 37, in map_expr_dag > result, = map_expr_dags(function, [expression], compress=compress) > File "/home/justin/Software/firedrake/src/ufl/ufl/corealg/map_dag.py", > line 84, in map_expr_dags > r = handlers[v._ufl_typecode_](v) > File > "/home/justin/Software/firedrake/src/ufl/ufl/algorithms/apply_restrictions.py", > line 102, in reference_value > g = self(f) > File > "/home/justin/Software/firedrake/src/ufl/ufl/corealg/multifunction.py", > line 95, in __call__ > return self._handlers[o._ufl_typecode_](o, *args) > File > "/home/justin/Software/firedrake/src/ufl/ufl/algorithms/apply_restrictions.py", > line 129, in coefficient > return self._require_restriction(o) > File > "/home/justin/Software/firedrake/src/ufl/ufl/algorithms/apply_restrictions.py", > line 59, in _require_restriction > "Discontinuous type %s must be restricted." % o._ufl_class_.__name__) > File "/home/justin/Software/firedrake/src/ufl/ufl/assertions.py", line > 51, in ufl_assert > error(*message) > File "/home/justin/Software/firedrake/src/ufl/ufl/log.py", line 167, in > error > raise self._exception_type(self._format_raw(*message)) > ufl.log.UFLException: Discontinuous type Coefficient must be restricted. > > I can't seem to figure out why. Attached is the code. Run as: > > python 2D_OS_boundary_ex1.py 100 50 0 0 0 > > Thanks, > Justin > > On Thu, Sep 29, 2016 at 4:49 AM, David Ham <David.Ham(a)imperial.ac.uk> > wrote: > > > > On Thu, 29 Sep 2016 at 10:34 Justin Chang <jychang48(a)gmail.com> wrote: > > Hi all, > > I have written here the Discontinuous Galerkin method for 1D > advection-diffusion equations using backward euler time discretization: > > mesh = UnitIntervalMesh(100) > V = FunctionSpace(mesh, 'DG', 1) > u = TrialFunction(V) > v = TestFunction(V) > > u0 = Function(V, name="initial") > sol = Function(V, name="solution") > > # Diffusivity > D = Constant(0.005) > > # Velocity > vel = as_vector((Constant(1.0),)) > > # Forcing function > f = interpolate(Expression(....),V) > > # Penalty > k = 1 > d = 1 > gamma = Constant((k+1)*(k+d)/d)*FacetArea(mesh)/avg(CellVolume(mesh)) # > based on the formula in > http://webpages.sdsmt.edu/~kshahbaz/shahbazi_05.pdf > > # Other parameters > n = FacetNormal(mesh) > h = CellSize(mesh) > vn = 0.5*(dot(vel,n) + abs(dot(vel,n))) > > # Weak form > a = v*u/dt*dx - dot(grad(v),vel*u)*dx + > dot(jump(v),vn('+')*u('+')-vn('-')*u('-'))*dS + dot(v, vn*u)*ds \ > + dot(grad(v),D*grad(u))*dx - dot(jump(v,n),avg(D*grad(u)))*dS - > dot(avg(D*grad(v)),jump(u,n))*dS \ > + gamma*dot(jump(v,n),jump(u,n))*dS > L = v*f*dx + v*u0/dt*dx > > > ... > > I have a couple questions: > > 1) If say I am solving for the transport of 4 chemical species, Normally I > would create a mixed function space like V_mixed = V*V*V*V, split u into > u1, u2, u3, and u4, write bilinear and linear forms for each of these > individual components and add them up. That is, a = a1 + a2 + a3 + a4, L = > L1 + L2 + L3 + L4. This can become quite a lot of code if I end up having > to solve ~20 chemical species eventually. The catch is that the velocity > and diffusivity will be the same for all species but the forcing function > and boundary conditions are going to be different. Is there a more > efficient of writing the bilinear and linear form instead of having to copy > and paste the a = ... and L = ... multiple times and changing the > trial/test functions? > > > a) If you are independently advecting lots of species (ie no reaction) > then you probably want to just solve the equations separately, rather than > making a massive mixed function space. > > b) You can avoid lots of retyping by just using normal Python functions to > construct forms. You could, for example, write a function which takes a > test and trial function, and maybe some BC information, and returns the > appropriate forms. > > > > 2) I see that in FEniCS's undocumented dg-advection-diffusion example, > they enforce the DirichletBC's strongly. Is it possible to enforce them > weakly like the SIPG or DG upwinding techniques for the pure diffusion or > advection equations? > > > Yes, just substitute boundary conditions for surface integrals in the > usual way. > > > > > Thanks, > Justin > > > > _______________________________________________ > firedrake mailing list > firedrake(a)imperial.ac.uk > https://mailman.ic.ac.uk/mailman/listinfo/firedrake > > > _______________________________________________ > firedrake mailing list > firedrake(a)imperial.ac.uk > https://mailman.ic.ac.uk/mailman/listinfo/firedrake > > > _______________________________________________ > firedrake mailing list > firedrake(a)imperial.ac.uk > https://mailman.ic.ac.uk/mailman/listinfo/firedrake > >

2 1

Re: [firedrake] DG for advection-diffusion
by Andrew McRae 05 Oct '16

05 Oct '16

Hi Justin, I don't think it's a good use of time for us to debug your full program (which runs an external python file at line 74, for example!) I suggest you strip down the program to its absolute core. If even this doesn't run then get back in touch with a minimal failing example. Thanks, Andrew On 5 October 2016 at 09:22, Justin Chang <jychang48(a)gmail.com> wrote: > Thanks, so the problem was that the diffusivity tensor was not inside the > avg() terms. That problem is fixed. > > However, when the code was modified as such, i get this error: > > Traceback (most recent call last): > File "2D_OS_boundary_ex1.py", line 313, in <module> > u0 = ADsolve(solver_D,tao_D,opt_D,u0) > File "2D_OS_boundary_ex1.py", line 292, in ADsolve > solver_O.solve(u0_O,b) > File "/home/justin/Software/firedrake/src/firedrake/firedrake/linear_solver.py", > line 141, in solve > raise ConvergenceError("LinearSolver failed to converge after %d > iterations with reason: %s", self.ksp.getIterationNumber(), > solving_utils.KSPReasons[r]) > firedrake.exceptions.ConvergenceError: ('LinearSolver failed to converge > after %d iterations with reason: %s', 0, 'DIVERGED_NANORINF') > > Is something wrong with my weak formulation? I even tried switching back > to weakly applied Boundaries but I still get the same error. > > Thanks, > Justin > > On Wed, Oct 5, 2016 at 3:05 AM, Andrew McRae <A.T.T.McRae(a)bath.ac.uk> > wrote: > >> Because all TrialFunctions, TestFunctions and Functions in interior facet >> integrals need to include restrictions (or use operators such as jump and >> avg) >> >> On 5 October 2016 at 08:32, Justin Chang <jychang48(a)gmail.com> wrote: >> >>> Okay thanks, >>> >>> I am attempting to do the diffusion part only at the moment under DG, >>> but I am getting these errors: >>> >>> $ python 2D_OS_boundary_ex1.py 100 50 0 0 0 -D_pc_type bjacobi >>> UFL:WARNING Discontinuous Lagrange element requested on quadrilateral, >>> creating DQ element. >>> UFL:WARNING Discontinuous Lagrange element requested on quadrilateral, >>> creating DQ element. >>> UFL:WARNING Discontinuous Lagrange element requested on quadrilateral, >>> creating DQ element. >>> UFL:ERROR Discontinuous type Coefficient must be restricted. >>> Traceback (most recent call last): >>> File "2D_OS_boundary_ex1.py", line 255, in <module> >>> A_D = assemble(a_D,bcs=bc_D) >>> File "/home/justin/Software/firedrake/src/firedrake/firedrake/assemble.py", >>> line 81, in assemble >>> inverse=inverse, mat_type=mat_type, appctx=appctx) >>> File "<decorator-gen-291>", line 2, in _assemble >>> File "/home/justin/Software/firedrake/src/firedrake/firedrake/utils.py", >>> line 62, in wrapper >>> return f(*args, **kwargs) >>> File "/home/justin/Software/firedrake/src/firedrake/firedrake/assemble.py", >>> line 120, in _assemble >>> inverse=inverse) >>> File "/home/justin/Software/firedrake/src/firedrake/firedrake/tsfc_interface.py", >>> line 189, in compile_form >>> number_map).kernels >>> File "/home/justin/Software/firedrake/src/PyOP2/pyop2/caching.py", >>> line 198, in __new__ >>> obj = make_obj() >>> File "/home/justin/Software/firedrake/src/PyOP2/pyop2/caching.py", >>> line 188, in make_obj >>> obj.__init__(*args, **kwargs) >>> File "/home/justin/Software/firedrake/src/firedrake/firedrake/tsfc_interface.py", >>> line 110, in __init__ >>> tree = tsfc_compile_form(form, prefix=name, parameters=parameters) >>> File "/home/justin/Software/firedrake/src/tsfc/tsfc/driver.py", line >>> 36, in compile_form >>> fd = ufl_utils.compute_form_data(form) >>> File "/home/justin/Software/firedrake/src/tsfc/tsfc/ufl_utils.py", >>> line 47, in compute_form_data >>> do_estimate_degrees=do_estimate_degrees, >>> File "/home/justin/Software/firedrake/src/ufl/ufl/algorithms/compute_form_data.py", >>> line 305, in compute_form_data >>> form = apply_restrictions(form) >>> File "/home/justin/Software/firedrake/src/ufl/ufl/algorithms/apply_restrictions.py", >>> line 211, in apply_restrictions >>> only_integral_type=integral_types) >>> File "/home/justin/Software/firedrake/src/ufl/ufl/algorithms/map_integrands.py", >>> line 58, in map_integrand_dags >>> form, only_integral_type) >>> File "/home/justin/Software/firedrake/src/ufl/ufl/algorithms/map_integrands.py", >>> line 39, in map_integrands >>> for itg in form.integrals()] >>> File "/home/justin/Software/firedrake/src/ufl/ufl/algorithms/map_integrands.py", >>> line 46, in map_integrands >>> return itg.reconstruct(function(itg.integrand())) >>> File "/home/justin/Software/firedrake/src/ufl/ufl/algorithms/map_integrands.py", >>> line 57, in <lambda> >>> return map_integrands(lambda expr: map_expr_dag(function, expr, >>> compress), >>> File "/home/justin/Software/firedrake/src/ufl/ufl/corealg/map_dag.py", >>> line 37, in map_expr_dag >>> result, = map_expr_dags(function, [expression], compress=compress) >>> File "/home/justin/Software/firedrake/src/ufl/ufl/corealg/map_dag.py", >>> line 84, in map_expr_dags >>> r = handlers[v._ufl_typecode_](v) >>> File "/home/justin/Software/firedrake/src/ufl/ufl/algorithms/apply_restrictions.py", >>> line 102, in reference_value >>> g = self(f) >>> File "/home/justin/Software/firedrake/src/ufl/ufl/corealg/multifunction.py", >>> line 95, in __call__ >>> return self._handlers[o._ufl_typecode_](o, *args) >>> File "/home/justin/Software/firedrake/src/ufl/ufl/algorithms/apply_restrictions.py", >>> line 129, in coefficient >>> return self._require_restriction(o) >>> File "/home/justin/Software/firedrake/src/ufl/ufl/algorithms/apply_restrictions.py", >>> line 59, in _require_restriction >>> "Discontinuous type %s must be restricted." % o._ufl_class_.__name__) >>> File "/home/justin/Software/firedrake/src/ufl/ufl/assertions.py", >>> line 51, in ufl_assert >>> error(*message) >>> File "/home/justin/Software/firedrake/src/ufl/ufl/log.py", line 167, >>> in error >>> raise self._exception_type(self._format_raw(*message)) >>> ufl.log.UFLException: Discontinuous type Coefficient must be restricted. >>> >>> I can't seem to figure out why. Attached is the code. Run as: >>> >>> python 2D_OS_boundary_ex1.py 100 50 0 0 0 >>> >>> Thanks, >>> Justin >>> >>> On Thu, Sep 29, 2016 at 4:49 AM, David Ham <David.Ham(a)imperial.ac.uk> >>> wrote: >>> >>>> >>>> >>>> On Thu, 29 Sep 2016 at 10:34 Justin Chang <jychang48(a)gmail.com> wrote: >>>> >>>>> Hi all, >>>>> >>>>> I have written here the Discontinuous Galerkin method for 1D >>>>> advection-diffusion equations using backward euler time discretization: >>>>> >>>>> mesh = UnitIntervalMesh(100) >>>>> V = FunctionSpace(mesh, 'DG', 1) >>>>> u = TrialFunction(V) >>>>> v = TestFunction(V) >>>>> >>>>> u0 = Function(V, name="initial") >>>>> sol = Function(V, name="solution") >>>>> >>>>> # Diffusivity >>>>> D = Constant(0.005) >>>>> >>>>> # Velocity >>>>> vel = as_vector((Constant(1.0),)) >>>>> >>>>> # Forcing function >>>>> f = interpolate(Expression(....),V) >>>>> >>>>> # Penalty >>>>> k = 1 >>>>> d = 1 >>>>> gamma = Constant((k+1)*(k+d)/d)*FacetArea(mesh)/avg(CellVolume(mesh)) >>>>> # based on the formula in http://webpages.sdsmt.edu/~ksh >>>>> ahbaz/shahbazi_05.pdf >>>>> >>>>> # Other parameters >>>>> n = FacetNormal(mesh) >>>>> h = CellSize(mesh) >>>>> vn = 0.5*(dot(vel,n) + abs(dot(vel,n))) >>>>> >>>>> # Weak form >>>>> a = v*u/dt*dx - dot(grad(v),vel*u)*dx + dot(jump(v),vn('+')*u('+')-vn('-')*u('-'))*dS >>>>> + dot(v, vn*u)*ds \ >>>>> + dot(grad(v),D*grad(u))*dx - dot(jump(v,n),avg(D*grad(u)))*dS - >>>>> dot(avg(D*grad(v)),jump(u,n))*dS \ >>>>> + gamma*dot(jump(v,n),jump(u,n))*dS >>>>> L = v*f*dx + v*u0/dt*dx >>>>> >>>>> >>>>> ... >>>>> >>>>> I have a couple questions: >>>>> >>>>> 1) If say I am solving for the transport of 4 chemical species, >>>>> Normally I would create a mixed function space like V_mixed = V*V*V*V, >>>>> split u into u1, u2, u3, and u4, write bilinear and linear forms for each >>>>> of these individual components and add them up. That is, a = a1 + a2 + a3 + >>>>> a4, L = L1 + L2 + L3 + L4. This can become quite a lot of code if I end up >>>>> having to solve ~20 chemical species eventually. The catch is that the >>>>> velocity and diffusivity will be the same for all species but the forcing >>>>> function and boundary conditions are going to be different. Is there a more >>>>> efficient of writing the bilinear and linear form instead of having to copy >>>>> and paste the a = ... and L = ... multiple times and changing the >>>>> trial/test functions? >>>>> >>>> >>>> a) If you are independently advecting lots of species (ie no reaction) >>>> then you probably want to just solve the equations separately, rather than >>>> making a massive mixed function space. >>>> >>>> b) You can avoid lots of retyping by just using normal Python functions >>>> to construct forms. You could, for example, write a function which takes a >>>> test and trial function, and maybe some BC information, and returns the >>>> appropriate forms. >>>> >>>> >>>>> >>>>> 2) I see that in FEniCS's undocumented dg-advection-diffusion example, >>>>> they enforce the DirichletBC's strongly. Is it possible to enforce them >>>>> weakly like the SIPG or DG upwinding techniques for the pure diffusion or >>>>> advection equations? >>>>> >>>> >>>> Yes, just substitute boundary conditions for surface integrals in the >>>> usual way. >>>> >>>> >>>> >>>>> >>>>> Thanks, >>>>> Justin >>>>> >>>>> >>>>> >>>> _______________________________________________ >>>> firedrake mailing list >>>> firedrake(a)imperial.ac.uk >>>> https://mailman.ic.ac.uk/mailman/listinfo/firedrake >>>> >>>> >>> >> >> _______________________________________________ >> firedrake mailing list >> firedrake(a)imperial.ac.uk >> https://mailman.ic.ac.uk/mailman/listinfo/firedrake >> >> >

2 1

Re: [firedrake] DG for advection-diffusion
by Andrew McRae 05 Oct '16

05 Oct '16

Because all TrialFunctions, TestFunctions and Functions in interior facet integrals need to include restrictions (or use operators such as jump and avg) On 5 October 2016 at 08:32, Justin Chang <jychang48(a)gmail.com> wrote: > Okay thanks, > > I am attempting to do the diffusion part only at the moment under DG, but > I am getting these errors: > > $ python 2D_OS_boundary_ex1.py 100 50 0 0 0 -D_pc_type bjacobi > UFL:WARNING Discontinuous Lagrange element requested on quadrilateral, > creating DQ element. > UFL:WARNING Discontinuous Lagrange element requested on quadrilateral, > creating DQ element. > UFL:WARNING Discontinuous Lagrange element requested on quadrilateral, > creating DQ element. > UFL:ERROR Discontinuous type Coefficient must be restricted. > Traceback (most recent call last): > File "2D_OS_boundary_ex1.py", line 255, in <module> > A_D = assemble(a_D,bcs=bc_D) > File "/home/justin/Software/firedrake/src/firedrake/firedrake/assemble.py", > line 81, in assemble > inverse=inverse, mat_type=mat_type, appctx=appctx) > File "<decorator-gen-291>", line 2, in _assemble > File "/home/justin/Software/firedrake/src/firedrake/firedrake/utils.py", > line 62, in wrapper > return f(*args, **kwargs) > File "/home/justin/Software/firedrake/src/firedrake/firedrake/assemble.py", > line 120, in _assemble > inverse=inverse) > File "/home/justin/Software/firedrake/src/firedrake/firedrake/tsfc_interface.py", > line 189, in compile_form > number_map).kernels > File "/home/justin/Software/firedrake/src/PyOP2/pyop2/caching.py", line > 198, in __new__ > obj = make_obj() > File "/home/justin/Software/firedrake/src/PyOP2/pyop2/caching.py", line > 188, in make_obj > obj.__init__(*args, **kwargs) > File "/home/justin/Software/firedrake/src/firedrake/firedrake/tsfc_interface.py", > line 110, in __init__ > tree = tsfc_compile_form(form, prefix=name, parameters=parameters) > File "/home/justin/Software/firedrake/src/tsfc/tsfc/driver.py", line > 36, in compile_form > fd = ufl_utils.compute_form_data(form) > File "/home/justin/Software/firedrake/src/tsfc/tsfc/ufl_utils.py", line > 47, in compute_form_data > do_estimate_degrees=do_estimate_degrees, > File "/home/justin/Software/firedrake/src/ufl/ufl/ > algorithms/compute_form_data.py", line 305, in compute_form_data > form = apply_restrictions(form) > File "/home/justin/Software/firedrake/src/ufl/ufl/ > algorithms/apply_restrictions.py", line 211, in apply_restrictions > only_integral_type=integral_types) > File "/home/justin/Software/firedrake/src/ufl/ufl/algorithms/map_integrands.py", > line 58, in map_integrand_dags > form, only_integral_type) > File "/home/justin/Software/firedrake/src/ufl/ufl/algorithms/map_integrands.py", > line 39, in map_integrands > for itg in form.integrals()] > File "/home/justin/Software/firedrake/src/ufl/ufl/algorithms/map_integrands.py", > line 46, in map_integrands > return itg.reconstruct(function(itg.integrand())) > File "/home/justin/Software/firedrake/src/ufl/ufl/algorithms/map_integrands.py", > line 57, in <lambda> > return map_integrands(lambda expr: map_expr_dag(function, expr, > compress), > File "/home/justin/Software/firedrake/src/ufl/ufl/corealg/map_dag.py", > line 37, in map_expr_dag > result, = map_expr_dags(function, [expression], compress=compress) > File "/home/justin/Software/firedrake/src/ufl/ufl/corealg/map_dag.py", > line 84, in map_expr_dags > r = handlers[v._ufl_typecode_](v) > File "/home/justin/Software/firedrake/src/ufl/ufl/ > algorithms/apply_restrictions.py", line 102, in reference_value > g = self(f) > File "/home/justin/Software/firedrake/src/ufl/ufl/corealg/multifunction.py", > line 95, in __call__ > return self._handlers[o._ufl_typecode_](o, *args) > File "/home/justin/Software/firedrake/src/ufl/ufl/ > algorithms/apply_restrictions.py", line 129, in coefficient > return self._require_restriction(o) > File "/home/justin/Software/firedrake/src/ufl/ufl/ > algorithms/apply_restrictions.py", line 59, in _require_restriction > "Discontinuous type %s must be restricted." % o._ufl_class_.__name__) > File "/home/justin/Software/firedrake/src/ufl/ufl/assertions.py", line > 51, in ufl_assert > error(*message) > File "/home/justin/Software/firedrake/src/ufl/ufl/log.py", line 167, in > error > raise self._exception_type(self._format_raw(*message)) > ufl.log.UFLException: Discontinuous type Coefficient must be restricted. > > I can't seem to figure out why. Attached is the code. Run as: > > python 2D_OS_boundary_ex1.py 100 50 0 0 0 > > Thanks, > Justin > > On Thu, Sep 29, 2016 at 4:49 AM, David Ham <David.Ham(a)imperial.ac.uk> > wrote: > >> >> >> On Thu, 29 Sep 2016 at 10:34 Justin Chang <jychang48(a)gmail.com> wrote: >> >>> Hi all, >>> >>> I have written here the Discontinuous Galerkin method for 1D >>> advection-diffusion equations using backward euler time discretization: >>> >>> mesh = UnitIntervalMesh(100) >>> V = FunctionSpace(mesh, 'DG', 1) >>> u = TrialFunction(V) >>> v = TestFunction(V) >>> >>> u0 = Function(V, name="initial") >>> sol = Function(V, name="solution") >>> >>> # Diffusivity >>> D = Constant(0.005) >>> >>> # Velocity >>> vel = as_vector((Constant(1.0),)) >>> >>> # Forcing function >>> f = interpolate(Expression(....),V) >>> >>> # Penalty >>> k = 1 >>> d = 1 >>> gamma = Constant((k+1)*(k+d)/d)*FacetArea(mesh)/avg(CellVolume(mesh)) # >>> based on the formula in http://webpages.sdsmt.edu/~ksh >>> ahbaz/shahbazi_05.pdf >>> >>> # Other parameters >>> n = FacetNormal(mesh) >>> h = CellSize(mesh) >>> vn = 0.5*(dot(vel,n) + abs(dot(vel,n))) >>> >>> # Weak form >>> a = v*u/dt*dx - dot(grad(v),vel*u)*dx + dot(jump(v),vn('+')*u('+')-vn('-')*u('-'))*dS >>> + dot(v, vn*u)*ds \ >>> + dot(grad(v),D*grad(u))*dx - dot(jump(v,n),avg(D*grad(u)))*dS - >>> dot(avg(D*grad(v)),jump(u,n))*dS \ >>> + gamma*dot(jump(v,n),jump(u,n))*dS >>> L = v*f*dx + v*u0/dt*dx >>> >>> >>> ... >>> >>> I have a couple questions: >>> >>> 1) If say I am solving for the transport of 4 chemical species, Normally >>> I would create a mixed function space like V_mixed = V*V*V*V, split u into >>> u1, u2, u3, and u4, write bilinear and linear forms for each of these >>> individual components and add them up. That is, a = a1 + a2 + a3 + a4, L = >>> L1 + L2 + L3 + L4. This can become quite a lot of code if I end up having >>> to solve ~20 chemical species eventually. The catch is that the velocity >>> and diffusivity will be the same for all species but the forcing function >>> and boundary conditions are going to be different. Is there a more >>> efficient of writing the bilinear and linear form instead of having to copy >>> and paste the a = ... and L = ... multiple times and changing the >>> trial/test functions? >>> >> >> a) If you are independently advecting lots of species (ie no reaction) >> then you probably want to just solve the equations separately, rather than >> making a massive mixed function space. >> >> b) You can avoid lots of retyping by just using normal Python functions >> to construct forms. You could, for example, write a function which takes a >> test and trial function, and maybe some BC information, and returns the >> appropriate forms. >> >> >>> >>> 2) I see that in FEniCS's undocumented dg-advection-diffusion example, >>> they enforce the DirichletBC's strongly. Is it possible to enforce them >>> weakly like the SIPG or DG upwinding techniques for the pure diffusion or >>> advection equations? >>> >> >> Yes, just substitute boundary conditions for surface integrals in the >> usual way. >> >> >> >>> >>> Thanks, >>> Justin >>> >>> >>> >> _______________________________________________ >> firedrake mailing list >> firedrake(a)imperial.ac.uk >> https://mailman.ic.ac.uk/mailman/listinfo/firedrake >> >> >

2 1

Re: [firedrake] DG for advection-diffusion
by David Ham 05 Oct '16

05 Oct '16

On Thu, 29 Sep 2016 at 10:34 Justin Chang <jychang48(a)gmail.com> wrote: > Hi all, > > I have written here the Discontinuous Galerkin method for 1D > advection-diffusion equations using backward euler time discretization: > > mesh = UnitIntervalMesh(100) > V = FunctionSpace(mesh, 'DG', 1) > u = TrialFunction(V) > v = TestFunction(V) > > u0 = Function(V, name="initial") > sol = Function(V, name="solution") > > # Diffusivity > D = Constant(0.005) > > # Velocity > vel = as_vector((Constant(1.0),)) > > # Forcing function > f = interpolate(Expression(....),V) > > # Penalty > k = 1 > d = 1 > gamma = Constant((k+1)*(k+d)/d)*FacetArea(mesh)/avg(CellVolume(mesh)) # > based on the formula in > http://webpages.sdsmt.edu/~kshahbaz/shahbazi_05.pdf > > # Other parameters > n = FacetNormal(mesh) > h = CellSize(mesh) > vn = 0.5*(dot(vel,n) + abs(dot(vel,n))) > > # Weak form > a = v*u/dt*dx - dot(grad(v),vel*u)*dx + > dot(jump(v),vn('+')*u('+')-vn('-')*u('-'))*dS + dot(v, vn*u)*ds \ > + dot(grad(v),D*grad(u))*dx - dot(jump(v,n),avg(D*grad(u)))*dS - > dot(avg(D*grad(v)),jump(u,n))*dS \ > + gamma*dot(jump(v,n),jump(u,n))*dS > L = v*f*dx + v*u0/dt*dx > > > ... > > I have a couple questions: > > 1) If say I am solving for the transport of 4 chemical species, Normally I > would create a mixed function space like V_mixed = V*V*V*V, split u into > u1, u2, u3, and u4, write bilinear and linear forms for each of these > individual components and add them up. That is, a = a1 + a2 + a3 + a4, L = > L1 + L2 + L3 + L4. This can become quite a lot of code if I end up having > to solve ~20 chemical species eventually. The catch is that the velocity > and diffusivity will be the same for all species but the forcing function > and boundary conditions are going to be different. Is there a more > efficient of writing the bilinear and linear form instead of having to copy > and paste the a = ... and L = ... multiple times and changing the > trial/test functions? > a) If you are independently advecting lots of species (ie no reaction) then you probably want to just solve the equations separately, rather than making a massive mixed function space. b) You can avoid lots of retyping by just using normal Python functions to construct forms. You could, for example, write a function which takes a test and trial function, and maybe some BC information, and returns the appropriate forms. > > 2) I see that in FEniCS's undocumented dg-advection-diffusion example, > they enforce the DirichletBC's strongly. Is it possible to enforce them > weakly like the SIPG or DG upwinding techniques for the pure diffusion or > advection equations? > Yes, just substitute boundary conditions for surface integrals in the usual way. > > Thanks, > Justin > > >

2 1

FSI coupling disagreement
by Tomasz Salwa [RPG] 29 Sep '16

29 Sep '16

Dear all (especially David and Lawrence), I work with a coupled fluid-structure interaction problem. I have a plain python code for verification. For a simple test setup (see attached file init_cond_py.pdf), I expect the same result from Firedrake as from my code. The result is very close, but slightly different. Grey plots show results from my code, coloured from Firedrake. The energy looks basically the same. However, in the inspection of convergence of energy with dt (expected second order convergence), Firedrake code results in a slight energy drift, which is unexpected. The uncoupled subsystems behave properly, so I suspect that the source of this discrepancy lies in the coupling terms. The relevant part is ... # given the definitions: V_f = FunctionSpace(mesh, "CG", 1) V_s = VectorFunctionSpace(mesh, "CG", 1) mixed_V = V_f * V_s trial_f, trial_s = TrialFunctions(mixed_V) v_f, v_s = TestFunctions(mixed_V) phi = Function(V_f) U = Function(V_s) n_normal = FacetNormal(mesh) # n - normal vector at the interface #I_s, I_f - indicator functions, 1(0) in solid, 0(1) in fluid, respectively n = (I_s("+") * I_f("-") - I_s("-") * I_f("+")) * n_normal("-") .... # the coupling terms are: ... + dot( avg(v_s), n ) * avg(trial_f) * dS ... + dot( avg(v_s), n ) * avg(phi) * dS ... - dot( n, avg(trial_s) ) * avg(v_f) * dS ... - dt * dot( n, avg(U) ) * avg(v_f) * dS Can it be that the energy drift is caused by the averages above or by a method Firedrake assembles those terms? I would welcome any suggestions, as I have run out of ideas what I may be doing wrong. The full code is here: https://bitbucket.org/mmtjs/lin_coupled_fd_2d with relevant part in lib/coupled_system.py Thank you, Tomasz

2 4