November 2015 - firedrake - imperial.ac.uk

Re: [firedrake] RT/BDM works but not TH/VMS
by Colin Cotter 13 Nov '15

13 Nov '15

Hi Justin, What is the form you are using for the VMS formulation? all the best --cjc On 13 November 2015 at 19:21, Justin Chang <jychang48(a)gmail.com> wrote: > Hi all, > > So given this Darcy problem: > > a = (dot(v, w) - p*div(w) - div(v)*q)*dx > L = f*q*dx > > where f = "12*pi*pi*sin(pi*x[0]*2)*sin(pi*x[1]*2)*sin(2*pi*x[2])" and with > these solver options: > > 'ksp_type': 'gmres', > #'ksp_monitor_true_residual': True, > 'pc_type': 'fieldsplit', > 'pc_fieldsplit_type': 'schur', > 'pc_fieldsplit_schur_fact_type': 'upper', > 'pc_fieldsplit_schur_precondition': 'selfp', > 'fieldsplit_0_ksp_type': 'preonly', > 'fieldsplit_0_pc_type': 'bjacobi', > 'fieldsplit_0_sub_pc_type': 'ilu', > 'fieldsplit_1_ksp_type': 'preonly', > 'fieldsplit_1_pc_type': 'hypre', > 'fieldsplit_1_pc_hypre_type': 'boomerang' > > And without specifying any boundaries (my forcing function f is chosen so > that the pressure is homogeneous on the boundary), RT0 and BDM works > beautifully. > > However, with Taylor-Hood elements and VMS (which is equal order CG1 plus > a least squares stabilization), the velocity works but the pressure > solution is screwed up. > > Why is that? > > Thanks, > Justin > -- http://www.imperial.ac.uk/people/colin.cotter www.cambridge.org/9781107663916

3 4

RT/BDM works but not TH/VMS
by Justin Chang 13 Nov '15

13 Nov '15

Hi all, So given this Darcy problem: a = (dot(v, w) - p*div(w) - div(v)*q)*dx L = f*q*dx where f = "12*pi*pi*sin(pi*x[0]*2)*sin(pi*x[1]*2)*sin(2*pi*x[2])" and with these solver options: 'ksp_type': 'gmres', #'ksp_monitor_true_residual': True, 'pc_type': 'fieldsplit', 'pc_fieldsplit_type': 'schur', 'pc_fieldsplit_schur_fact_type': 'upper', 'pc_fieldsplit_schur_precondition': 'selfp', 'fieldsplit_0_ksp_type': 'preonly', 'fieldsplit_0_pc_type': 'bjacobi', 'fieldsplit_0_sub_pc_type': 'ilu', 'fieldsplit_1_ksp_type': 'preonly', 'fieldsplit_1_pc_type': 'hypre', 'fieldsplit_1_pc_hypre_type': 'boomerang' And without specifying any boundaries (my forcing function f is chosen so that the pressure is homogeneous on the boundary), RT0 and BDM works beautifully. However, with Taylor-Hood elements and VMS (which is equal order CG1 plus a least squares stabilization), the velocity works but the pressure solution is screwed up. Why is that? Thanks, Justin

1 0

Re: [firedrake] Projecting CG1 function into DG0 function
by David Ham 13 Nov '15

13 Nov '15

On Fri, 13 Nov 2015 at 16:29 Justin Chang <jychang48(a)gmail.com> wrote: > Yeah I found out about errornorm recently. Have been using that now. > > But as for the original issue... would D.assign(assemble(div(q) * e * dx) > suffice? > I'm afraid not. assemble(div(q)*e*dx) is \int_\Omega div(q) e dx \forall e \in E where E is the basis for DG0. This is not the same thing as the projection of div(q) into DG0. In fact, for the very specific case of DG0, this will differ from the correct answer in each element by a factor of the volume of that element. This is why dividing by \int_\Omega e dx \forall e \in E is necessary to achieve the correct answer. To go into slightly more detail, note that the e_i \in E is: / 1 if x \in T_i (the i-th element) e_i(x) = | \ 0 otherwise. Hence: \int_\Omega \div(q) e_i dx == \int_{T_i} div(q)) dx and \int_\Omega e_i dx == \int_{T_i} dx == vol(T_i) but it must be the case that: \int_T_i D_i e_i dx == \int_{T_i} div(q)) dx and note: \int_T_i D_i e_i dx == D_i \int_T_i dx (since D_i is just a scalar, not a function of x) hence by rearranging: D_i = \int_\Omega \div(q) e_i dx/ \int_\Omega e_i dx and hence: D = \int_\Omega \div(q) e_i dx/ \int_\Omega e_i dx \forall e_i \in E > > On Fri, Nov 13, 2015 at 8:27 AM, David Ham <David.Ham(a)imperial.ac.uk> > wrote: > >> >> >> On Thu, 12 Nov 2015 at 21:48 Justin Chang <jychang48(a)gmail.com> wrote: >> >>> Hi guys, >>> >>> So I was attempting to do something similar: >>> >>> L2_error_norm = norm(assemble(v*(u-u_exact)*dx)/assemble(v*dx)) >>> >>> where v is test function, u is FE solution, and u_exact is the >>> analytical solution. All of which are CG1 space. I get an error saying >>> "fl.log.UFLException: Division by non-scalar is undefined" >>> >>> Know what's up? >>> >> >> Yes, that's either a missing feature or a bug depending on how you look >> at it. However I think you are just confused in what you are trying to do. >> >> What you probably meant was: >> >> errornorm(u, u_exact) >> >> (errornorm is a Firedrake-provided function). >> >> What errornorm is basically doing is: >> >> assemble((u-u_exact)**2*dx))**0.5 >> >> Note that this is a zero-form (no test or trial function; the answer is a >> scalar) so there are no solves, or divisions by lumped mass. >> >> errornorm actually does somewhat more sophisticated things such as >> increasing the polynomial order of the functions in order to avoid certain >> sorts of error, so it's typically better to use errornorm than to do it >> yourself. >> >> Regards, >> >> David >> >> >> >>> >>> Thanks, >>> Justin >>> >>> On Thu, Nov 12, 2015 at 4:33 AM, Lawrence Mitchell < >>> lawrence.mitchell(a)imperial.ac.uk> wrote: >>> >>>> On 12/11/15 11:27, Justin Chang wrote: >>>> > David, >>>> > >>>> > So if D.assign(assemble(div(q) * e * dx)/assemble(e * dx)) returns >>>> > cell-wise div(q), what's the denominator "/assemble(e * dx)" for? >>>> >>>> >>>> It's just the normal FE L2 projection: >>>> >>>> you want: >>>> >>>> u = div(q) >>>> >>>> So you hit both sides with a test function and integrate: >>>> >>>> u*e*dx = div(q)*e*dx >>>> >>>> Where u is a trial function in DG0, e is a test function in DG and q >>>> is in whereever. >>>> >>>> But, as David points out, the DG0 mass-matrix is completely diagonal. >>>> The values on the diagonal are just obtained by assembling the lone >>>> test function: >>>> >>>> diag = assemble(e*dx) >>>> >>>> But now, because the matrix is diagonal, you can solve the linear >>>> system pointwise by doing a pointwise division: >>>> >>>> rhs = assemble(div(q)*e*dx) >>>> >>>> solution = assemble(rhs/diag) >>>> >>>> >>>> Cheers, >>>> >>>> Lawrence >>>> >>>> >>>> _______________________________________________ >>>> 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 >> >> >

1 0

Re: [firedrake] Projecting CG1 function into DG0 function
by David Ham 13 Nov '15

13 Nov '15

On Thu, 12 Nov 2015 at 21:48 Justin Chang <jychang48(a)gmail.com> wrote: > Hi guys, > > So I was attempting to do something similar: > > L2_error_norm = norm(assemble(v*(u-u_exact)*dx)/assemble(v*dx)) > > where v is test function, u is FE solution, and u_exact is the analytical > solution. All of which are CG1 space. I get an error saying > "fl.log.UFLException: Division by non-scalar is undefined" > > Know what's up? > Yes, that's either a missing feature or a bug depending on how you look at it. However I think you are just confused in what you are trying to do. What you probably meant was: errornorm(u, u_exact) (errornorm is a Firedrake-provided function). What errornorm is basically doing is: assemble((u-u_exact)**2*dx))**0.5 Note that this is a zero-form (no test or trial function; the answer is a scalar) so there are no solves, or divisions by lumped mass. errornorm actually does somewhat more sophisticated things such as increasing the polynomial order of the functions in order to avoid certain sorts of error, so it's typically better to use errornorm than to do it yourself. Regards, David > > Thanks, > Justin > > On Thu, Nov 12, 2015 at 4:33 AM, Lawrence Mitchell < > lawrence.mitchell(a)imperial.ac.uk> wrote: > >> On 12/11/15 11:27, Justin Chang wrote: >> > David, >> > >> > So if D.assign(assemble(div(q) * e * dx)/assemble(e * dx)) returns >> > cell-wise div(q), what's the denominator "/assemble(e * dx)" for? >> >> >> It's just the normal FE L2 projection: >> >> you want: >> >> u = div(q) >> >> So you hit both sides with a test function and integrate: >> >> u*e*dx = div(q)*e*dx >> >> Where u is a trial function in DG0, e is a test function in DG and q >> is in whereever. >> >> But, as David points out, the DG0 mass-matrix is completely diagonal. >> The values on the diagonal are just obtained by assembling the lone >> test function: >> >> diag = assemble(e*dx) >> >> But now, because the matrix is diagonal, you can solve the linear >> system pointwise by doing a pointwise division: >> >> rhs = assemble(div(q)*e*dx) >> >> solution = assemble(rhs/diag) >> >> >> Cheers, >> >> Lawrence >> >> >> _______________________________________________ >> firedrake mailing list >> firedrake(a)imperial.ac.uk >> https://mailman.ic.ac.uk/mailman/listinfo/firedrake >> >> >

2 1

Output for time-stepping?
by Justin Chang 13 Nov '15

13 Nov '15

Hi all, So I am running simulations which have many time steps and normally I would save the files as .pvd via outfile << File(...). This creates a new file for each time level However, if I need to solve a large geochemistry system of equations (~20 unknowns) over 100+ time steps, I would prefer not to have 2000+ pvd files sitting my file system/laptop. Is there a way to store the solution of all time levels into a single file? Thanks, Justin

4 7

Re: [firedrake] Problem with Jacobian
by David Ham 13 Nov '15

13 Nov '15

On Mon, 9 Nov 2015 at 08:31 Buesing, Henrik <HBuesing(a)eonerc.rwth-aachen.de> wrote: > > If you change the definition of derivative in firedrake/ufl_expr.py to > > > take an optional coefficient_derivatives argument: > > > > > > def derivative(F, u, du=None, coefficient_derivatives=None): > > > ... > > > > > > > > > And just pass that value through to the ufl.derivative call: > > > > > > ufl.derivative(..., coefficient_derivatives=coefficient_derivatives) > > > > > > If you do this in the firedrake source directly, you will be able to > > > commit your change in the local git repository. You can send us these > > > changes by forking firedrake on github and proposing a pull-request that > > > way. > > Hmm... Okay I did this. Then I did a "make" in the firedrake directory. > But then I'm still getting a > > > > *TypeError: derivative() takes at most 3 arguments (4 given)* > > > > So how can I test this? Thank you! > > > > It depends a little on how you originally installed firedrake. If you > installed firedrake by running firedrake-install (i.e. not with > --developer) then you change to the firedrake/src directory and (with the > virtualenv active) type: > > > > pip install --upgrade firedrake/ > > > > If you have installed firedrake some other way, please say how and we can > tell you how to build. > > > > *[Buesing, Henrik] I followed the „Installing from individual components” > manual. I cloned the firedrake repository* > > git clone https://github.com/firedrakeproject/firedrake.git > > *Then I did make in the firedrake directory.* > > cd firedrake > > export PYTHONPATH=$PWD:$PYTHONPATH > > make > > In this setup you shouldn't have even needed to type "make". My initial suspicion is that this is an actual issue with the new code. Can you please push your changes to a fork on github so we can see what you've got there? Regards, David > > > *Henrik* >

3 3

Matrix multiplication in bilinear form
by Anna Kalogirou 13 Nov '15

13 Nov '15

Dear all, Can I use an assembled matrix into a bilinear form? E.g. I want to write something like mu1 = Function(V) mu = TrialFunction(V) v = TestFunction(V) Q = assemble(v*dx) lhs = inner(grad(mu),grad(v))*dx + Q*mu*dx rhs = .... solve(lhs, mu1, rhs) Shall/can I just assemble each term in the lhs independently, and then add/multiply the resulting matrices? Best, Anna.

2 14

Re: [firedrake] New Firedrake dependency - libspatialindex
by David Ham 13 Nov '15

13 Nov '15

Thanks for catching this, Justin. Miklos has made a fix to the install script accordingly. Cheers, David On Thu, 29 Oct 2015 at 19:33 Justin Chang <jychang48(a)gmail.com> wrote: > Hi Lawrence, > > Yeah the issue was that /usr/local/include was not shown originally. I ran > code-select --install and everything runs fine now. > > Thanks, > Justin > > On Thu, Oct 29, 2015 at 3:21 AM, Lawrence Mitchell < > lawrence.mitchell(a)imperial.ac.uk> wrote: > >> -----BEGIN PGP SIGNED MESSAGE----- >> Hash: SHA1 >> >> On 29/10/15 07:45, Justin Chang wrote: >> > I just did a fresh install on my Macbook (also upgraded to El >> > Capitan) but I am still getting errors. Attached is the log file. >> > >> > Can you guys tell what's going on? >> >> >> firedrake/spatialindex.cpp:259:10: fatal error: >> 'spatialindex/SpatialIndex.h' file not found >> #include "spatialindex/SpatialIndex.h" >> >> Suggests that although spatialindex appears to be installed, your >> compiler does not find it. This header file should be in >> /usr/local/include >> >> So what do the following say? >> >> $ ls /usr/local/include/spatialindex >> $ mpicc -x c -v -E /dev/null >> >> The latter should print out a section: >> >> #include <...> search starts here: >> LIST OF DIRECTORIES >> End of search list >> >> If /usr/local/include doesn't appear in this list of directories can >> you try: >> >> $ xcode-select --install >> >> And then check again. >> >> Cheers, >> >> Lawrence >> -----BEGIN PGP SIGNATURE----- >> Version: GnuPG v2.0.22 (GNU/Linux) >> >> iQEcBAEBAgAGBQJWMeU2AAoJECOc1kQ8PEYv6lUH/R1JnxVnmIDheeqQ4Z0xWuSq >> ePznzKtGb4cQPRddBxksC46miENCdeXF3zhOEAfqGV5eaAum/7xW7FgzbiPHy+dl >> 9V1RyKnzQm+WXHS0dfWuqa0QOJITnJivzyMZxm7Rt5Ayqa/S5AK2yq2fsI399JI8 >> Nwt53/yFjRvy/pja7YhCTMfb/gquz82ZP/EvtU5RWcHEh+Fr1MQ///GzJVRDGgGk >> ldSxdawDJW8sAdmUkdtLmwEj1tfRXRVtrdKg0fzgFWjEzhFbzOH8lUNavoKsLJHu >> SpIjFvc1BEbR67hzsJugeRN9ZWNJhwLn1X5Hu1CIhYfM+IR0kx9SnkJKMYHjc3I= >> =0S8+ >> -----END PGP SIGNATURE----- >> >> _______________________________________________ >> firedrake mailing list >> firedrake(a)imperial.ac.uk >> https://mailman.ic.ac.uk/mailman/listinfo/firedrake >> > >

3 3

Re: [firedrake] Projecting CG1 function into DG0 function
by Andrew McRae 12 Nov '15

12 Nov '15

Is 'norm' some symbolic UFL thing? On 12 November 2015 at 21:48, Justin Chang <jychang48(a)gmail.com> wrote: > Hi guys, > > So I was attempting to do something similar: > > L2_error_norm = norm(assemble(v*(u-u_exact)*dx)/assemble(v*dx)) > > where v is test function, u is FE solution, and u_exact is the analytical > solution. All of which are CG1 space. I get an error saying > "fl.log.UFLException: Division by non-scalar is undefined" > > Know what's up? > > Thanks, > Justin > > On Thu, Nov 12, 2015 at 4:33 AM, Lawrence Mitchell < > lawrence.mitchell(a)imperial.ac.uk> wrote: > >> On 12/11/15 11:27, Justin Chang wrote: >> > David, >> > >> > So if D.assign(assemble(div(q) * e * dx)/assemble(e * dx)) returns >> > cell-wise div(q), what's the denominator "/assemble(e * dx)" for? >> >> >> It's just the normal FE L2 projection: >> >> you want: >> >> u = div(q) >> >> So you hit both sides with a test function and integrate: >> >> u*e*dx = div(q)*e*dx >> >> Where u is a trial function in DG0, e is a test function in DG and q >> is in whereever. >> >> But, as David points out, the DG0 mass-matrix is completely diagonal. >> The values on the diagonal are just obtained by assembling the lone >> test function: >> >> diag = assemble(e*dx) >> >> But now, because the matrix is diagonal, you can solve the linear >> system pointwise by doing a pointwise division: >> >> rhs = assemble(div(q)*e*dx) >> >> solution = assemble(rhs/diag) >> >> >> Cheers, >> >> Lawrence >> >> >> _______________________________________________ >> firedrake mailing list >> firedrake(a)imperial.ac.uk >> https://mailman.ic.ac.uk/mailman/listinfo/firedrake >> >> >

2 1

Re: [firedrake] Projecting CG1 function into DG0 function
by David Ham 12 Nov '15

12 Nov '15

Basically: d = TrialFunction(D.function_space()) e = TestFunction(D.function_space()) solve(d * e * dx == div(q) * e * dx, D) One might have to think a little about what value this gives at the boundary, but I think it's OK. On Thu, 12 Nov 2015 at 11:00 Justin Chang <jychang48(a)gmail.com> wrote: > Hi all, > > Perhaps this may be a simple question, but say I have this bilinear and > linear form: > > a = grad(u)*D*grad(v)*dx > L = F*v*dx > > where u,v is trial/test function on CG1 space, and D and F are > coefficients that live in DG0 space. > > Say I have a velocity vector function q (CG1) on the same mesh and want D > to be the element-wise divergence of said velocity. > > How would I formulate D in Firedrake? > > Thanks, > Justin >

3 4