October 2016 - firedrake - imperial.ac.uk

Re: [firedrake] MeshHierarchy problem
by Colin Cotter 29 Oct '16

29 Oct '16

Hi Wei, How critical is it to have periodic bcs for your problem? We could alternately pick slip bcs which avoid this problem. All the best Cjc On Saturday, 29 October 2016, Pan, Wei <wei.pan(a)imperial.ac.uk> wrote: > Hi Lawrence, > > > Many thanks for the reply. > > > > You could maybe do it (a bit of hack) by building a TorusMesh, then > creating a mesh hierarchy of that, and then using the code in > firedrake/utility_meshes.py to do create your periodic meshes and put them > back into the hierarchy (ugh). > > > I do need to compare solutions for different mesh refinements. So you mean > I should create the TorusMesh hierarchy and build the PeriodicSquareMesh > from the torus meshes? I'm unsure what you mean by putting the periodic > meshes back into the hierarchy. How do I do that? Is there a way to create > a 1-to-1 mapping between a torus mesh and periodic square mesh? > > > I will submit a github issue. > > > Kind regards > > Wei > ------------------------------ > *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 Lawrence Mitchell <lawrence.mitchell(a)imperial.ac.uk > <javascript:_e(%7B%7D,'cvml','lawrence.mitchell(a)imperial.ac.uk');>> > *Sent:* 28 October 2016 14:15:52 > *To:* firedrake > *Subject:* Re: [firedrake] MeshHierarchy problem > > > > On 28 Oct 2016, at 13:24, Pan, Wei <wei.pan(a)imperial.ac.uk > <javascript:_e(%7B%7D,'cvml','wei.pan(a)imperial.ac.uk');>> wrote: > > > > Hi, > > > > I'm currently using a mixed formulation for solving the Euler equation > with periodic boundary conditions. So right now I am using a > PeriodicSquareMesh. I also need different refinements of the same mesh. I > was made aware of the MeshHierarchy function. > > Do you need to transfer solutions between the different refinements? If > not, just make meshes of the appropriate resolution. Otherwise, read on... > > > However MeshHierarchy does not seem to work with PeriodicSquareMesh, I > get an error message saying > > > > "ValueError: total size of new array must be unchanged." > > > I think this is because we assume when building meshes from a DM object > that the coordinate field lives in a piecewise linear space. However, this > is not the case for periodic meshes. We rely on PETSc refinement to refine > the mesh, but for periodic meshes, PETSc doesn't have a full view of the > mesh (it only really knows the topology). > > You could maybe do it (a bit of hack) by building a TorusMesh, then > creating a mesh hierarchy of that, and then using the code in > firedrake/utility_meshes.py to do create your periodic meshes and put them > back into the hierarchy (ugh). > > In any case, can you please report a github issue describing what you > tried to do and what didn't work? > > Cheers, > > Lawrence > _______________________________________________ > firedrake mailing list > firedrake(a)imperial.ac.uk > <javascript:_e(%7B%7D,'cvml','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

1 0

MeshHierarchy problem
by Pan, Wei 29 Oct '16

29 Oct '16

Hi, I'm currently using a mixed formulation for solving the Euler equation with periodic boundary conditions. So right now I am using a PeriodicSquareMesh. I also need different refinements of the same mesh. I was made aware of the MeshHierarchy function. However MeshHierarchy does not seem to work with PeriodicSquareMesh, I get an error message saying "ValueError: total size of new array must be unchanged." Is it possible to do periodic boundary conditions on SquareMesh? Many thanks Wei

2 2

function node subset
by Tomasz Salwa [RPG] 27 Oct '16

27 Oct '16

Dear all, How do I properly construct the node subset given function values? I'd like to get something like this: a = Function(V) # ... initialize a with some values ... b = Function(V) b.assign(1/a, subset=["where a != 0"] ) Thank you, Tomasz

1 0

mapping between two meshes
by Tomasz Salwa [RPG] 26 Oct '16

26 Oct '16

Dear all, I would like to check the convergence of my solution as I refine the mesh. Say, I have two solutions in CG1 space, f_coarse on mesh_coarse and f_fine on mesh_fine, which is a refined mesh_coarse. Can I map/interpolate/project f_fine onto mesh_coarse in Firedrake? Thank you, Tomasz

3 5

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

24 Oct '16

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? cheers --cjc On 7 October 2016 at 14:39, Anna Kalogirou <A.Kalogirou(a)leeds.ac.uk> wrote: > Hi Lawrence, > > I don’t use any preconditioning in either of the two sets of code. > > Anna. > > > On 7 Oct 2016, at 14:33, Lawrence Mitchell <lawrence.mitchell@imperial. > ac.uk> wrote: > > > Hi Anna, > > > >> On 6 Oct 2016, at 18:42, Anna Kalogirou <A.Kalogirou(a)leeds.ac.uk> > wrote: > >> > >> Sorry for the delay, I was away from my desk all day today. I attach > the two log files. > >> > >> Best, Anna. > >> > >> <log_file_mixed.txt><log_file.txt> > > > > OK great. We can see that both options have quite a large room for > improvement in the linear solver. > > > > The original code does ~6000 KSPSolve calls, but calls PCApply ~395000 > times. This tells me you're taking around 65 iterations per solve. Which > is probably quite a lot unless your problem is very stiff. > > > > The new code (mixed system) is worse, you do only ~2000 KSPSolve calls, > but call PCApply almost ~3000000 times, or around ~1500 iterations per > solve. This is terrible :). > > > > What preconditioning, if any, are you using for these solves (in both > cases?). As Colin says, we need to figure out what the linear system looks > like to choose an appropriate preconditioner. > > > > 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 > -- http://www.imperial.ac.uk/people/colin.cotter www.cambridge.org/9781107663916

4 5

visualising 2nd order element
by Schwedes, Tobias 21 Oct '16

21 Oct '16

Hi all, I would like to visualize the result of simulations based on 2nd order elements. However, using plot(u) only gives a linear interpolation. I know this from fenics and the workaround is to interpolate u onto a finer mesh and plot the outcome then. Does that work in firedrake, and if so how do I do it? Thanks a lot in advance for your help! All the best, Tobias

2 5

Files don't exist
by Justin Chang 20 Oct '16

20 Oct '16

Hi all, I am looking through the "Support for matrix-free operator application and preconditioning" section of the online documentation but none of the links to the demo python files work. I know it's extremely minor, but are they typos or simply non-existent? Thanks, Justin

1 0

assign different values on subdomains
by Schwedes, Tobias 19 Oct '16

19 Oct '16

Hi all, I would like to assign different values of some function u on subdomains of a 2-D domain. More precisely, I want to have u equal to 1 for y > x and equal to 0 otherwise. For that reason I guess one needs to be able to define subdomains. Besides the Issue on https://github.com/firedrakeproject/firedrake/issues/723 I didn't find any helpful information online. Can someone help out please? I attached the fenics version of what I want, hoping that might help! All the best, Tobias

3 3

Re: [firedrake] assign different values on subdomains
by Colin Cotter 18 Oct '16

18 Oct '16

Hi Tobias, This can alternatively be done by making u to be a discontinuous field and projecting your expression. On Tuesday, 18 October 2016, Schwedes, Tobias <t.schwedes14(a)imperial.ac.uk> wrote: > Hi all, > > > I would like to assign different values of some function u on subdomains > of a 2-D domain. More precisely, I want to have u equal to 1 for y > x > and equal to 0 otherwise. For that reason I guess one needs to be able to > define subdomains. Besides the Issue on > > > https://github.com/firedrakeproject/firedrake/issues/723 > > > I didn't find any helpful information online. Can someone help out please? > I attached the fenics version of what I want, hoping that might help! > > > All the best, > > Tobias > -- http://www.imperial.ac.uk/people/colin.cotter www.cambridge.org/9781107663916

1 0

Re: [firedrake] Python expression for loading files
by Andrew McRae 13 Oct '16

13 Oct '16

Just to revisit this one last time, the main reason this was slow is because you were accessing x_DG.dat.data_ro[j,0] and x_DG.dat.data_ro[j,1] inside the innermost, i, loop (twice each!). This came with *a lot* of Python overhead. Since these quantities are independent of i, you could have defined temporary variables outside the i loop: x0 = x_DG.dat.data_ro[j,0] x1 = x_DG.dat.data_ro[j,1] and used these variables within the loop. Just doing this makes the code ~4x faster (75% saving) than what was in your email. Getting rid of the innermost i loop and using array-sized operations, as I did before, then gives the rest of the benefit (an extra factor of 20x or more...) On 12 October 2016 at 21:18, Justin Chang <jychang48(a)gmail.com> wrote: > Yes that works brilliantly. Thanks Andrew! > > On Wed, Oct 12, 2016 at 2:37 PM, Andrew McRae <A.T.T.McRae(a)bath.ac.uk> > wrote: > >> FWIW, this already gives a *massive* speedup: >> >> for j in range(usize): >> num = np.zeros(600) >> den = np.zeros(600) >> >> num[:] = datas[:]/(4*np.power(dataxs[:] - x_DG.dat.data_ro[j, 0], 2) + >> np.power(datays[:] - x_DG.dat.data_ro[j, 1], 2)) >> >> den[:] = 1/(4*np.power(dataxs[:] - x_DG.dat.data_ro[j, 0], 2) + >> np.power(datays[:] - x_DG.dat.data_ro[j, 1], 2)) >> >> u.dat.data[j] = 1e-12*np.exp(sum(num)/sum(den)) >> >> There might still be a smarter way to do it. >> >> On 12 October 2016 at 20:23, Andrew McRae <A.T.T.McRae(a)bath.ac.uk> wrote: >> >>> Justin, >>> >>> Your loops are in standard Python. This is likely to be slow. >>> >>> To give a simple case, if you want to calculate a = b + c, where these >>> are 1000x1000 NumPy arrays, your current code looks like >>> >>> for j in range(1000): >>> for i in range(1000): >>> a[i, j] = b[i, j] + c[i, j] >>> >>> David is suggesting that you replace this with the "NumPy operation" >>> >>> a[:, :] = b[:, :] + c[:, :] >>> >>> which will be much faster. In this simple case you could just write a = >>> b + c, but in your real case you'll probably need to use slice notation. >>> >>> >>> >>> >>> On 12 October 2016 at 19:44, Justin Chang <jychang48(a)gmail.com> wrote: >>> >>>> David, >>>> >>>> I tried to follow your steps, but I ended up with something that's even >>>> worse performance wise. Clearly I am not doing something right. >>>> >>>> Attached is my code. Am I not utilizing numpy operations correctly? >>>> >>>> Thanks, >>>> Justin >>>> >>>> On Wed, Oct 12, 2016 at 5:13 AM, David Ham <David.Ham(a)imperial.ac.uk> >>>> wrote: >>>> >>>>> Hi Justin, >>>>> >>>>> Yes. >>>>> >>>>> 1. Make a DG0 "coordinate" field by projecting the coordinates into >>>>> VectorFunctionSpace(mesh, DG, "0") call it, say X_DG >>>>> 2. Use numpy operations to collectively interpolate your input data to >>>>> the points given by X_DG.dat.data_ro >>>>> 3. Set u.dat.data[:] to those interpolated values. >>>>> >>>>> This will result in going through the python stack O(1) times instead >>>>> of O(meshsize) which should be waaay faster. >>>>> >>>>> On Wed, 12 Oct 2016 at 11:00 Justin Chang <jychang48(a)gmail.com> wrote: >>>>> >>>>>> Miklos I am not sure I follow what you're saying. >>>>>> >>>>>> Also, my code is pretty slow lol, is there a "more efficient" way of >>>>>> doing this or is this something I will just have to deal with? >>>>>> >>>>>> On Wed, Oct 12, 2016 at 4:46 AM, Homolya, Miklós < >>>>>> m.homolya14(a)imperial.ac.uk> wrote: >>>>>> >>>>>> Maybe you need to call the __init__ is the superclass? >>>>>> ------------------------------ >>>>>> *From:* firedrake-bounces(a)imperial.ac.uk < >>>>>> firedrake-bounces(a)imperial.ac.uk> on behalf of Justin Chang < >>>>>> jychang48(a)gmail.com> >>>>>> *Sent:* 12 October 2016 10:45:11 >>>>>> *To:* firedrake >>>>>> *Subject:* Re: [firedrake] Python expression for loading files >>>>>> >>>>>> Nevermind, I fixed it by adding this line under __init__ >>>>>> >>>>>> self._user_args = [] >>>>>> >>>>>> Code works now. Thanks again for the help >>>>>> >>>>>> Justin >>>>>> >>>>>> On Wed, Oct 12, 2016 at 4:38 AM, Justin Chang <jychang48(a)gmail.com> >>>>>> wrote: >>>>>> >>>>>> I tried that, but I now get this error: >>>>>> >>>>>> Traceback (most recent call last): >>>>>> File "readData.py", line 25, in <module> >>>>>> u.interpolate(createData()) >>>>>> File "/home/justin/Software/firedrake/src/firedrake/firedrake/function.py", >>>>>> line 328, in interpolate >>>>>> return interpolation.interpolate(expression, self, subset=subset) >>>>>> File "/home/justin/Software/firedrake/src/firedrake/firedrake/interpolation.py", >>>>>> line 36, in interpolate >>>>>> return Interpolator(expr, V, subset=subset).interpolate() >>>>>> File "/home/justin/Software/firedrake/src/firedrake/firedrake/interpolation.py", >>>>>> line 56, in __init__ >>>>>> self.callable = make_interpolator(expr, V, subset) >>>>>> File "/home/justin/Software/firedrake/src/firedrake/firedrake/interpolation.py", >>>>>> line 116, in make_interpolator >>>>>> loops.append(_interpolator(V, f.dat, expr, subset)) >>>>>> File "/home/justin/Software/firedrake/src/firedrake/firedrake/interpolation.py", >>>>>> line 166, in _interpolator >>>>>> kernel, oriented, coefficients = compile_python_kernel(expr, >>>>>> to_pts, to_element, V, coords) >>>>>> File "/home/justin/Software/firedrake/src/firedrake/firedrake/interpolation.py", >>>>>> line 310, in compile_python_kernel >>>>>> for _, arg in expression._user_args: >>>>>> AttributeError: 'createData' object has no attribute '_user_args' >>>>>> >>>>>> >>>>>> Here's my code and text files >>>>>> >>>>>> Thanks, >>>>>> Justin >>>>>> >>>>>> On Wed, Oct 12, 2016 at 4:36 AM, Homolya, Miklós < >>>>>> m.homolya14(a)imperial.ac.uk> wrote: >>>>>> >>>>>> def value_shape(self): >>>>>> return (1,) >>>>>> >>>>>> Change (1,) to () >>>>>> >>>>>> >>>>>> >>>>>> On Wed, Oct 12, 2016 at 10:35 AM +0100, "Justin Chang" < >>>>>> jychang48(a)gmail.com> wrote: >>>>>> >>>>>> How do I make the expression return a scalar? >>>>>> >>>>>> On Wed, Oct 12, 2016 at 4:28 AM, Homolya, Miklós < >>>>>> m.homolya14(a)imperial.ac.uk> wrote: >>>>>> >>>>>> Your function space Q is scalar (), while your expression is a 1D >>>>>> vector (1,). Change one of them to the other. >>>>>> >>>>>> >>>>>> >>>>>> On Wed, Oct 12, 2016 at 10:05 AM +0100, "Justin Chang" < >>>>>> jychang48(a)gmail.com> wrote: >>>>>> >>>>>> Dear all, >>>>>> >>>>>> I am attempting to create an Expression which uses data from some >>>>>> text files. This is my code: >>>>>> >>>>>> from firedrake import * >>>>>> import numpy as np >>>>>> >>>>>> mesh = UnitSquareMesh(100,100) >>>>>> Q = FunctionSpace(mesh,'DG',0) >>>>>> >>>>>> file_prefix = 'perm600' >>>>>> class createData(Expression): >>>>>> def __init__(self): >>>>>> self.s = np.loadtxt(file_prefix+'_s') >>>>>> self.xs = np.loadtxt(file_prefix+'_xs') >>>>>> self.ys = np.loadtxt(file_prefix+'_ys') >>>>>> def eval(self, val, x): >>>>>> numerator = 0 >>>>>> denominator = 0 >>>>>> for i in range(0,600): >>>>>> numerator += self.s[i]/((self.xs[i]-x[0])** >>>>>> 2+(self.ys[i]-x[1])**2) >>>>>> denominator += 1/((self.xs[i]-x[0])**2+(self.ys[i]-x[1])**2) >>>>>> ssum = numerator/denominator >>>>>> val[0] = 1e-12*exp(ssum) >>>>>> def value_shape(self): >>>>>> return (1,) >>>>>> >>>>>> u = Function(Q) >>>>>> u.interpolate(createData()) >>>>>> >>>>>> >>>>>> However, I get this error: >>>>>> >>>>>> Traceback (most recent call last): >>>>>> File "readData.py", line 25, in <module> >>>>>> u.interpolate(createData()) >>>>>> File "/home/justin/Software/firedrake/src/firedrake/firedrake/function.py", >>>>>> line 328, in interpolate >>>>>> return interpolation.interpolate(expression, self, subset=subset) >>>>>> File "/home/justin/Software/firedrake/src/firedrake/firedrake/interpolation.py", >>>>>> line 36, in interpolate >>>>>> return Interpolator(expr, V, subset=subset).interpolate() >>>>>> File "/home/justin/Software/firedrake/src/firedrake/firedrake/interpolation.py", >>>>>> line 56, in __init__ >>>>>> self.callable = make_interpolator(expr, V, subset) >>>>>> File "/home/justin/Software/firedrake/src/firedrake/firedrake/interpolation.py", >>>>>> line 116, in make_interpolator >>>>>> loops.append(_interpolator(V, f.dat, expr, subset)) >>>>>> File "/home/justin/Software/firedrake/src/firedrake/firedrake/interpolation.py", >>>>>> line 152, in _interpolator >>>>>> % (len(expr.ufl_shape), len(V.ufl_element().value_shape()))) >>>>>> RuntimeError: Rank mismatch: Expression rank 1, FunctionSpace rank 0 >>>>>> >>>>>> What's going on here? I know you guys generally recommend using UFL >>>>>> for interpolating expressions, but I am not sure how I would implement the >>>>>> above in a UFL way since I have to read data from a textfile. >>>>>> >>>>>> Any help appreciated :) >>>>>> >>>>>> 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 >> >> >

1 0