June 2016 - firedrake - imperial.ac.uk

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

25 Jun '16

Onno did attach rectmesh.msh... Onno, is there a reason why you didn't use one of the built-in meshes? cheers --cjc On 25 June 2016 at 16:18, Francis Poulin <fpoulin(a)uwaterloo.ca> wrote: > Hello Onno, > > I just tried running your code but see that I can't because I don't have > the rectmesh.msh. Might you be able to send that so we can see exactly > what you see? > > Francis > > P.S. I would also be curious to see the gmsh file that generated the msh > file if possible. > > > > ------------------ > Francis Poulin > Associate Professor > Department of Applied Mathematics > University of Waterloo > > email: fpoulin(a)uwaterloo.ca > Web: https://uwaterloo.ca/poulin-research-group/ > Telephone: +1 519 888 4567 x32637 > > ------------------------------ > *From:* firedrake-bounces(a)imperial.ac.uk [firedrake-bounces(a)imperial.ac.uk] > on behalf of Onno Bokhove [O.Bokhove(a)leeds.ac.uk] > *Sent:* Saturday, June 25, 2016 11:05 AM > *To:* firedrake > *Subject:* [firedrake] upwind example remark; extension to CN in time > > Hi FDs, > > > (i) The upwind example is odd in that the real boundaries are treated > > in a special way rather than setting the outside value in > > the same numerical flux as anywhere else; > > my question is how this can be fixed? > > So no special terms at the in- and outflow, just specify > > the outside velocity used. This is imo the proper DG way not was is done > > now. > > > (ii) Attached my attempt to do a time dependent advection > > with Crank Nicolson, which works but does not yet > > give the right answer, in part due to (i). > > > [I don't understand why I have to make > > my theta variable that a constant but other wise it did not work.] > > Any suggestions, corrections, also on other silly things done by me? > > Next step is to extend to system with 3 equations; > hence the few extra but unused variables. > > Thanks, from midnight Seoul in old style Korean guesthouse. > > Onno > ------------------------------ > -- http://www.imperial.ac.uk/people/colin.cotter www.cambridge.org/9781107663916

2 1

Re: [firedrake] missing 'comm'
by Colin Cotter 25 Jun '16

25 Jun '16

OK, now it is not finding some libraries: In file included from /usr/include/spatialindex/SpatialIndex.h:30:0, from /usr/include/spatialindex/capi/sidx_config.h:56, from /usr/include/spatialindex/capi/sidx_api.h:34, from firedrake/spatialindex.c:245: /usr/include/spatialindex/tools/Tools.h:59:20: fatal error: iostream: No such file or directory #include <iostream> ^ compilation terminated. error: command 'mpicc' failed with exit status 1 all the best --cjc On 25 June 2016 at 15:41, Homolya, Miklós <m.homolya14(a)imperial.ac.uk> wrote: > I have seen such problems already, so it might be systematic, but I'm not > sure what are the exact circumstances that trigger this problem. > > > Try: > > > make clean all > > > in the Firedrake source directory. > ------------------------------ > *From:* firedrake-bounces(a)imperial.ac.uk <firedrake-bounces(a)imperial.ac.uk> > on behalf of Colin Cotter <colin.cotter(a)imperial.ac.uk> > *Sent:* 25 June 2016 14:54:41 > *To:* firedrake > *Subject:* [firedrake] missing 'comm' > > Dear Firedrakers, > Upon calling > > mesh = IcosahedralSphereMesh(1.0,3,2) > > I get an error, > > Traceback (most recent call last): > File "dcmip_3-1_linear_advection_triangle_bdfm1.py", line 31, in <module> > extrusion_type="radial") > File "<decorator-gen-274>", line 2, in ExtrudedMesh > File "/home/cjc1/firedrake-dev/firedrake/src/PyOP2/pyop2/profiling.py", > line 59, in wrapper > return f(*args, **kwargs) > File > "/home/cjc1/firedrake-dev/firedrake/src/firedrake/firedrake/mesh.py", line > 1232, in ExtrudedMesh > mesh.init() > File > "/home/cjc1/firedrake-dev/firedrake/src/firedrake/firedrake/mesh.py", line > 831, in init > self._callback(self) > File > "/home/cjc1/firedrake-dev/firedrake/src/firedrake/firedrake/mesh.py", line > 1165, in callback > self.topology.init() > File "<decorator-gen-272>", line 2, in init > File "/home/cjc1/firedrake-dev/firedrake/src/PyOP2/pyop2/profiling.py", > line 59, in wrapper > return f(*args, **kwargs) > File > "/home/cjc1/firedrake-dev/firedrake/src/firedrake/firedrake/mesh.py", line > 407, in init > self._callback(self) > File > "/home/cjc1/firedrake-dev/firedrake/src/firedrake/firedrake/mesh.py", line > 375, in callback > dmplex.mark_entity_classes(self._plex) > File "firedrake/dmplex.pyx", line 632, in > firedrake.dmplex.mark_entity_classes (firedrake/dmplex.c:9423) > AttributeError: 'module' object has no attribute 'comm' > > Is this another petsc4py issue? > > cheers > --cjc > > -- > http://www.imperial.ac.uk/people/colin.cotter > > www.cambridge.org/9781107663916 > > > -- http://www.imperial.ac.uk/people/colin.cotter www.cambridge.org/9781107663916

3 2

missing 'comm'
by Colin Cotter 25 Jun '16

25 Jun '16

Dear Firedrakers, Upon calling mesh = IcosahedralSphereMesh(1.0,3,2) I get an error, Traceback (most recent call last): File "dcmip_3-1_linear_advection_triangle_bdfm1.py", line 31, in <module> extrusion_type="radial") File "<decorator-gen-274>", line 2, in ExtrudedMesh File "/home/cjc1/firedrake-dev/firedrake/src/PyOP2/pyop2/profiling.py", line 59, in wrapper return f(*args, **kwargs) File "/home/cjc1/firedrake-dev/firedrake/src/firedrake/firedrake/mesh.py", line 1232, in ExtrudedMesh mesh.init() File "/home/cjc1/firedrake-dev/firedrake/src/firedrake/firedrake/mesh.py", line 831, in init self._callback(self) File "/home/cjc1/firedrake-dev/firedrake/src/firedrake/firedrake/mesh.py", line 1165, in callback self.topology.init() File "<decorator-gen-272>", line 2, in init File "/home/cjc1/firedrake-dev/firedrake/src/PyOP2/pyop2/profiling.py", line 59, in wrapper return f(*args, **kwargs) File "/home/cjc1/firedrake-dev/firedrake/src/firedrake/firedrake/mesh.py", line 407, in init self._callback(self) File "/home/cjc1/firedrake-dev/firedrake/src/firedrake/firedrake/mesh.py", line 375, in callback dmplex.mark_entity_classes(self._plex) File "firedrake/dmplex.pyx", line 632, in firedrake.dmplex.mark_entity_classes (firedrake/dmplex.c:9423) AttributeError: 'module' object has no attribute 'comm' Is this another petsc4py issue? cheers --cjc -- http://www.imperial.ac.uk/people/colin.cotter www.cambridge.org/9781107663916

2 1

install script problem when building spatialindex
by Colin Cotter 25 Jun '16

25 Jun '16

Dear FIredrakers, I currently get this message from the install/update script: Cloning libspatialindex Enter passphrase for key '/home/cjc1/.ssh/id_rsa': Successfully cloned libspatialindex branch master. HEAD is now at 34fd22b Merge pull request #82 from miklos1/master glibtoolize or libtoolize not found. Giving up! Traceback (most recent call last): File "firedrake-install", line 903, in <module> build_and_install_libspatialindex() File "firedrake-install", line 617, in build_and_install_libspatialindex check_call(["./autogen.sh"]) File "firedrake-install", line 211, in check_call subprocess.check_call(arguments, env=env) File "/usr/lib/python2.7/subprocess.py", line 540, in check_call raise CalledProcessError(retcode, cmd) subprocess.CalledProcessError: Command '['./autogen.sh']' returned non-zero exit status 1 What can I do to fix this? cheers --cjc -- http://www.imperial.ac.uk/people/colin.cotter www.cambridge.org/9781107663916

2 2

instability of a vortex in firedrakeQG
by Francis Poulin 17 Jun '16

17 Jun '16

Hello, I recently posted a video on youtube that shows the evolution of the perturbation Potential Vorticity for a Gaussian vortex. Initially you see a tripole structure and then a mode 2 structure that is well predicted by the linear theory. This is all done using the QG firedrake code that we've been developing. https://youtu.be/6YLPT7glvlU One thing that is less than ideal is that I am using white noise for the initial perturbations. Someone mentioned that I should probably solve the Poisson problem where white noise is the RHS, and that should clean things up. I will try and do that but it seems that early on, even after the perturbations establish themselves you can still see the grid. In this case I am using a 400x400 rectangular grid. Do people think that an unstructured mesh might make the grid less apparent? Maybe part of the problem is with how it's being plotted in VisIt, I'm not sure but I thought I would ask people's advice. Cheers, Francis ------------------ Francis Poulin Associate Professor Department of Applied Mathematics University of Waterloo email: fpoulin(a)uwaterloo.ca Web: https://uwaterloo.ca/poulin-research-group/ Telephone: +1 519 888 4567 x32637

2 3

Re: [firedrake] removing nan's?
by David Ham 17 Jun '16

17 Jun '16

Or: qb.dat.data[:] = np.nan_to_num(qb.dat.data_ro) On Fri, 17 Jun 2016 at 08:47 Andrew McRae <A.T.T.McRae(a)bath.ac.uk> wrote: > I think the more direct approach is to edit the numpy array qb.dat.data. > > E.g. something like qb.dat.data[np.isnan(qb.dat.data)] = 0 > > On 17 June 2016 at 03:52, Francis Poulin <fpoulin(a)uwaterloo.ca> wrote: > >> Hello, >> >> I want to define an initial condition for a field that I am pretty sure >> will have nan's. It does in python so I'm pretty sure that it will in >> firedrake as well. >> >> In python I know how to remove that nan's by doing the following. >> >> qb[np.isnan(qb)] = 0. >> >> Is there a way of doing this in firedrake? >> >> Cheers, Francis >> >> >> ------------------ >> Francis Poulin >> Associate Professor >> Department of Applied Mathematics >> University of Waterloo >> >> email: fpoulin(a)uwaterloo.ca >> Web: https://uwaterloo.ca/poulin-research-group/ >> Telephone: +1 519 888 4567 x32637 >> >> >

3 2

Re: [firedrake] feet wet
by Colin Cotter 17 Jun '16

17 Jun '16

Hi Onno, If you check out the QG demo, it uses an RK scheme with DG for pv. I've been intending to supply a demo for DG transient advection fora while, maybe over the summer I'll find time. Cheers cjc On Friday, 17 June 2016, Onno Bokhove <O.Bokhove(a)leeds.ac.uk> wrote: > Hi FD's > > > Trying to get my feet wet. Simple query. > > > Trying to make a DG conservation law/system code, > > starting with linear advection, backward Euler in steps. > > > How do I make the code snippet below forward Euler > > or Crank-Nicolson? I used all the old time step terms in L > > but that does not work. (Now commented out -thet factor terms) > > Also, why does the nonlinear solver not work? > > I think it is because the printing can't be done in DG0 > > so there needs to assign between the DG solution of the > > prime unknown. > > > Would be nice to have a standard DG example for > > U_t + div F = 0 > > forward Euler or RK3. > > > Thanks for the reply as then I can keep playing in the plane! > > > Onno > > > P.S. Code below > > > tend = 4.211 > > timestep = tend/40 > > > V = FunctionSpace(tmesh, "DG", 0) > > W = VectorFunctionSpace(tmesh, "CG", 1) > > > # Advection velocity > > velocity = Expression(("1.0", "0.0")) > > u0 = project(velocity, W) > > #Inflow BC > > inflow = Expression("(x[0] < 0.0) ? 1.0 : 0") > > inflow2 = Expression("(x[0] < 0.5) ? 1.0 : 0") > > inflow3 = Expression("(x[0] < 1.0) ? 1.0 : 0") > > #Inflow Variable > > D0 = Function(V) > > D0.interpolate(inflow) > > m1 = Function(V) > > m1.interpolate(inflow2) > > m2 = Function(V) > > m2.interpolate(inflow3) > > > #Initial Condition > > ic = Expression("(x[0] < .5) ? 1.0 : 0") > > ic2 = Expression("(x[0] < .5) ? 1.0 : 0") > > ic3 = Expression("(x[0] < .5) ? 1.0 : 0") > > > #Set up element normal > > n = FacetNormal(tmesh) > > # > > un = 0.5*(dot(u0, n) + abs(dot(u0, n))) > > > # Previous timestep for BE timestep > > D_ = Function(V) > > m1_ = Function(V) > > m2_ = Function(V) > > > D_.interpolate(ic) > > m1_.interpolate(ic2) > > m2_.interpolate(ic3) > > > D = TrialFunction(V) > > phi = TestFunction(V) > > # > > # Setting up Bilinear form > > # > > thet = 0; > > a1 = ((D/timestep)*phi - (1-thet)*D*dot(u0, grad(phi)))*dx # > (1-thet) interior term at new time step > > a2 = (1-thet)*dot(jump(phi), un('+')*D('+') - un('-')*D('-'))*dS # > Internal flux > > a3 = (1-thet)*dot(phi, un*D)*(ds(112)) # outflow > at top wall > > a = a1 + a2 + a3 > > # > > L1 = (D_/timestep)*phi*dx - D0*phi*dot(u0, n)*(ds(114)) # > previous timestep & given inflow at bottom wall > > # L2 = thet*(D_*dot(u0, grad(phi)))*dx # thet > times interior term at old time stepL > > # L3 = -thet*dot(jump(phi), un('+')*D_('+') - un('-')*D_('-'))*dS # > Internal flux > > # L4 = -thet*dot(phi, un*D_)*(ds(112)) # outflow > at top wall > > L = L1 > > F = a-L1 > > > out = Function(V) > > outfile = File("./Results/swe2d.pvd") > > # Write IC to file > > # outfile << D_ > > D_.assign(out) > > outfile.write(out) > > # outfile.wrote(project(out, W, name="D")) > > t = 0.0 > > while (t <= tend): > > #Solve problem > > solve(a == L, out) > > #solve(F == 0, out) > > #Update previous timestep > > D_.assign(out) > > #Update time > > t+=timestep > > # outfile << out > > outfile.write(out) > > # outfile.write(project(out, W, name="D")) > >> >> -- http://www.imperial.ac.uk/people/colin.cotter www.cambridge.org/9781107663916

1 0

Re: [firedrake] removing nan's?
by Andrew McRae 17 Jun '16

17 Jun '16

I think the more direct approach is to edit the numpy array qb.dat.data. E.g. something like qb.dat.data[np.isnan(qb.dat.data)] = 0 On 17 June 2016 at 03:52, Francis Poulin <fpoulin(a)uwaterloo.ca> wrote: > Hello, > > I want to define an initial condition for a field that I am pretty sure > will have nan's. It does in python so I'm pretty sure that it will in > firedrake as well. > > In python I know how to remove that nan's by doing the following. > > qb[np.isnan(qb)] = 0. > > Is there a way of doing this in firedrake? > > Cheers, Francis > > > ------------------ > Francis Poulin > Associate Professor > Department of Applied Mathematics > University of Waterloo > > email: fpoulin(a)uwaterloo.ca > Web: https://uwaterloo.ca/poulin-research-group/ > Telephone: +1 519 888 4567 x32637 > >

1 0

Re: [firedrake] removing nan's?
by Colin Cotter 17 Jun '16

17 Jun '16

Hi Francis, If you can predict the situations under which you get divide by zero, then you can use the ?: operator in an Expression string or Condition in ufl. All the best Cjc On Friday, 17 June 2016, Francis Poulin <fpoulin(a)uwaterloo.ca> wrote: > Hello, > > I want to define an initial condition for a field that I am pretty sure > will have nan's. It does in python so I'm pretty sure that it will in > firedrake as well. > > In python I know how to remove that nan's by doing the following. > > qb[np.isnan(qb)] = 0. > > Is there a way of doing this in firedrake? > > Cheers, Francis > > > ------------------ > Francis Poulin > Associate Professor > Department of Applied Mathematics > University of Waterloo > > email: fpoulin(a)uwaterloo.ca > <javascript:_e(%7B%7D,'cvml','fpoulin(a)uwaterloo.ca');> > Web: https://uwaterloo.ca/poulin-research-group/ > Telephone: +1 519 888 4567 x32637 > > -- http://www.imperial.ac.uk/people/colin.cotter www.cambridge.org/9781107663916

1 0

removing nan's?
by Francis Poulin 17 Jun '16

17 Jun '16

Hello, I want to define an initial condition for a field that I am pretty sure will have nan's. It does in python so I'm pretty sure that it will in firedrake as well. In python I know how to remove that nan's by doing the following. qb[np.isnan(qb)] = 0. Is there a way of doing this in firedrake? Cheers, Francis ------------------ Francis Poulin Associate Professor Department of Applied Mathematics University of Waterloo email: fpoulin(a)uwaterloo.ca Web: https://uwaterloo.ca/poulin-research-group/ Telephone: +1 519 888 4567 x32637

1 0