December 2015 - firedrake - imperial.ac.uk

PETSc via petsc4py/PyOP2 really slow
by Justin Chang 04 Jan '16

04 Jan '16

Hi all, If I had a purely firedrake code, everything runs fast on my Macbook. But when I need to make direct calls to the TaoSolver via petsc4py or PETSc calls via PyOP2, the program is extremely slow. This was not the case as of several commits ago (like probably two months). I suspect it has something to do with not adding any optimization flags to the PETSc configure options? The PETSc developers suggest adding more than simply --with-debugging=0. For instance, if using GNU compilers one could use COPTFLAGS='-O3 -march=native' It seems to me that whatever compiler optimization flags were used for Firedrake/PyOP2 function calls are not used for direct usage of PETSc function calls? Thanks, Justin

2 2

Applying bcs when using LinearSolver
by Justin Chang 04 Jan '16

04 Jan '16

Hi all, Say I have this: a = dot(grad(u),grad(v))*dx L = f*v*dx bcs = DirichletBC(...) A = assemble(a) b = assemble(L) solver = LinearSolver(A, ...) solver.solve(x,b) If I want to apply the DirichletBC's, I would use "A = assemble(a,bcs=bcs)." Is it necessary to also do this for b? Would it have any affect? Thanks, Justin

3 4

Re: [firedrake] Animation
by Andrew McRae 22 Dec '15

22 Dec '15

Yes, but outfile is only defined once, and gets written to every 10th step. I removed the output mess from the code I just sent... On 22 December 2015 at 14:04, Colin Cotter <colin.cotter(a)imperial.ac.uk> wrote: > Hi Angwenyi, > You're still defining "output" inside the time loop. > > all the best > --cjc > > On 22 December 2015 at 13:18, Angwenyi David <kipkoej(a)gmail.com> wrote: > >> >> On 22 Dec 2015, at 14:12, Andrew McRae <A.T.T.McRae(a)bath.ac.uk> wrote: >> >> Send your complete code please. >> >> >> >>> from firedrake import * >> >> >>> nx = 10 >> >> >>> ny = 10 >> >> >>> nz = 5 >> >> >>> mesh = BoxMesh(nx, ny, nz, 2.0, 2.0, 1.0) >> >> >>> V = FunctionSpace(mesh, 'Lagrange', 1) >> >> >>> p = Function(V, name="p") >> >> >>> phi = Function(V, name="phi") >> >> >>> >> >> >>> u = TrialFunction(V) >> >> >>> v = TestFunction(V) >> >> >>> outfile = File("out.pvd") >> >> >>> outfile << phi >> >> >>> bcval = Constant(0.0) >> >> >>> bc = DirichletBC(V, bcval, 1) >> >> >>> T = 10. >> >> >>> dt = 0.001 >> >> >>> t = 0 >> >> >>> step = 0 >> >> >>> lump_mass = True >> >> >>> while t <= T: >> >> ... step += 1 >> >> ... bcval.assign(sin(2*pi*5*t)) >> >> ... phi -= dt / 2 * p >> >> ... if lump_mass: >> >> ... p += assemble(dt * inner(nabla_grad(v), nabla_grad(phi))*dx) >> / assemble(v*dx) >> >> ... bc.apply(p) >> >> ... else: >> >> ... solve(u * v * dx == v * p * dx + dt * inner(grad(v), >> grad(phi)) * dx,p, bcs=bc, solver_parameters={'ksp_type': 'cg','pc_type': >> 'sor','pc_sor_symmetric': True}) >> >> ... phi -= dt / 2 * p >> >> ... t += dt >> >> ... if step % 10 == 0: >> >> ... outfile << phi >> >> ... File("waveEqn.pvd")<< phi >> >> ... output = File("waveEqn.pvd") >> >> ... output << phi >> >> ... del output >> >> ... >> >> >> >> >> On 22 December 2015 at 13:10, Angwenyi David <kipkoej(a)gmail.com> wrote: >> >>> Dear Cotter, dear Andrew, >>> >>> I get one .vtu file in both cases. >>> >>> Angwenyi >>> >>> On 22 Dec 2015, at 13:11, Colin Cotter <colin.cotter(a)imperial.ac.uk> >>> wrote: >>> >>> Dear Angwenyi, >>> The output should be one .pvd file, but many .vtu files, one per dump. >>> Then, you open the .pvd file in Paraview, set up your visualisation and >>> press the play button. >>> >>> all the best >>> --cjc >>> >>> On 22 December 2015 at 12:06, Angwenyi David <kipkoej(a)gmail.com> wrote: >>> >>>> Dear Cotter, >>>> >>>> I did so, but got a single file. I thought I could get as many files as >>>> timesteps; which files I could open all at once in Paraview to generate a >>>> video. Any direction? >>>> >>>> Angwenyi >>>> >>>> On 22 Dec 2015, at 12:10, Colin Cotter <colin.cotter(a)imperial.ac.uk> >>>> wrote: >>>> >>>> del output >>>> >>>> >>>> >>> >>> >>> -- >>> 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 >> >> >> > > > -- > http://www.imperial.ac.uk/people/colin.cotter > > www.cambridge.org/9781107663916 > > >

2 1

Re: [firedrake] Animation
by Colin Cotter 22 Dec '15

22 Dec '15

Hi Angwenyi, You're still defining "output" inside the time loop. all the best --cjc On 22 December 2015 at 13:18, Angwenyi David <kipkoej(a)gmail.com> wrote: > > On 22 Dec 2015, at 14:12, Andrew McRae <A.T.T.McRae(a)bath.ac.uk> wrote: > > Send your complete code please. > > > >>> from firedrake import * > > >>> nx = 10 > > >>> ny = 10 > > >>> nz = 5 > > >>> mesh = BoxMesh(nx, ny, nz, 2.0, 2.0, 1.0) > > >>> V = FunctionSpace(mesh, 'Lagrange', 1) > > >>> p = Function(V, name="p") > > >>> phi = Function(V, name="phi") > > >>> > > >>> u = TrialFunction(V) > > >>> v = TestFunction(V) > > >>> outfile = File("out.pvd") > > >>> outfile << phi > > >>> bcval = Constant(0.0) > > >>> bc = DirichletBC(V, bcval, 1) > > >>> T = 10. > > >>> dt = 0.001 > > >>> t = 0 > > >>> step = 0 > > >>> lump_mass = True > > >>> while t <= T: > > ... step += 1 > > ... bcval.assign(sin(2*pi*5*t)) > > ... phi -= dt / 2 * p > > ... if lump_mass: > > ... p += assemble(dt * inner(nabla_grad(v), nabla_grad(phi))*dx) / > assemble(v*dx) > > ... bc.apply(p) > > ... else: > > ... solve(u * v * dx == v * p * dx + dt * inner(grad(v), > grad(phi)) * dx,p, bcs=bc, solver_parameters={'ksp_type': 'cg','pc_type': > 'sor','pc_sor_symmetric': True}) > > ... phi -= dt / 2 * p > > ... t += dt > > ... if step % 10 == 0: > > ... outfile << phi > > ... File("waveEqn.pvd")<< phi > > ... output = File("waveEqn.pvd") > > ... output << phi > > ... del output > > ... > > > > > On 22 December 2015 at 13:10, Angwenyi David <kipkoej(a)gmail.com> wrote: > >> Dear Cotter, dear Andrew, >> >> I get one .vtu file in both cases. >> >> Angwenyi >> >> On 22 Dec 2015, at 13:11, Colin Cotter <colin.cotter(a)imperial.ac.uk> >> wrote: >> >> Dear Angwenyi, >> The output should be one .pvd file, but many .vtu files, one per dump. >> Then, you open the .pvd file in Paraview, set up your visualisation and >> press the play button. >> >> all the best >> --cjc >> >> On 22 December 2015 at 12:06, Angwenyi David <kipkoej(a)gmail.com> wrote: >> >>> Dear Cotter, >>> >>> I did so, but got a single file. I thought I could get as many files as >>> timesteps; which files I could open all at once in Paraview to generate a >>> video. Any direction? >>> >>> Angwenyi >>> >>> On 22 Dec 2015, at 12:10, Colin Cotter <colin.cotter(a)imperial.ac.uk> >>> wrote: >>> >>> del output >>> >>> >>> >> >> >> -- >> 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 > > > -- http://www.imperial.ac.uk/people/colin.cotter www.cambridge.org/9781107663916

1 0

Re: [firedrake] Animation
by Andrew McRae 22 Dec '15

22 Dec '15

I believe there should be multiple out... '.vtu' files for your code. If not, please try the attached code. Best, Andrew On 22 December 2015 at 13:18, Angwenyi David <kipkoej(a)gmail.com> wrote: > > On 22 Dec 2015, at 14:12, Andrew McRae <A.T.T.McRae(a)bath.ac.uk> wrote: > > Send your complete code please. > > > >>> from firedrake import * > > >>> nx = 10 > > >>> ny = 10 > > >>> nz = 5 > > >>> mesh = BoxMesh(nx, ny, nz, 2.0, 2.0, 1.0) > > >>> V = FunctionSpace(mesh, 'Lagrange', 1) > > >>> p = Function(V, name="p") > > >>> phi = Function(V, name="phi") > > >>> > > >>> u = TrialFunction(V) > > >>> v = TestFunction(V) > > >>> outfile = File("out.pvd") > > >>> outfile << phi > > >>> bcval = Constant(0.0) > > >>> bc = DirichletBC(V, bcval, 1) > > >>> T = 10. > > >>> dt = 0.001 > > >>> t = 0 > > >>> step = 0 > > >>> lump_mass = True > > >>> while t <= T: > > ... step += 1 > > ... bcval.assign(sin(2*pi*5*t)) > > ... phi -= dt / 2 * p > > ... if lump_mass: > > ... p += assemble(dt * inner(nabla_grad(v), nabla_grad(phi))*dx) / > assemble(v*dx) > > ... bc.apply(p) > > ... else: > > ... solve(u * v * dx == v * p * dx + dt * inner(grad(v), > grad(phi)) * dx,p, bcs=bc, solver_parameters={'ksp_type': 'cg','pc_type': > 'sor','pc_sor_symmetric': True}) > > ... phi -= dt / 2 * p > > ... t += dt > > ... if step % 10 == 0: > > ... outfile << phi > > ... File("waveEqn.pvd")<< phi > > ... output = File("waveEqn.pvd") > > ... output << phi > > ... del output > > ... > > > > > On 22 December 2015 at 13:10, Angwenyi David <kipkoej(a)gmail.com> wrote: > >> Dear Cotter, dear Andrew, >> >> I get one .vtu file in both cases. >> >> Angwenyi >> >> On 22 Dec 2015, at 13:11, Colin Cotter <colin.cotter(a)imperial.ac.uk> >> wrote: >> >> Dear Angwenyi, >> The output should be one .pvd file, but many .vtu files, one per dump. >> Then, you open the .pvd file in Paraview, set up your visualisation and >> press the play button. >> >> all the best >> --cjc >> >> On 22 December 2015 at 12:06, Angwenyi David <kipkoej(a)gmail.com> wrote: >> >>> Dear Cotter, >>> >>> I did so, but got a single file. I thought I could get as many files as >>> timesteps; which files I could open all at once in Paraview to generate a >>> video. Any direction? >>> >>> Angwenyi >>> >>> On 22 Dec 2015, at 12:10, Colin Cotter <colin.cotter(a)imperial.ac.uk> >>> wrote: >>> >>> del output >>> >>> >>> >> >> >> -- >> 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 > > >

1 0

Re: [firedrake] Animation
by Andrew McRae 22 Dec '15

22 Dec '15

Send your complete code please. On 22 December 2015 at 13:10, Angwenyi David <kipkoej(a)gmail.com> wrote: > Dear Cotter, dear Andrew, > > I get one .vtu file in both cases. > > Angwenyi > > On 22 Dec 2015, at 13:11, Colin Cotter <colin.cotter(a)imperial.ac.uk> > wrote: > > Dear Angwenyi, > The output should be one .pvd file, but many .vtu files, one per dump. > Then, you open the .pvd file in Paraview, set up your visualisation and > press the play button. > > all the best > --cjc > > On 22 December 2015 at 12:06, Angwenyi David <kipkoej(a)gmail.com> wrote: > >> Dear Cotter, >> >> I did so, but got a single file. I thought I could get as many files as >> timesteps; which files I could open all at once in Paraview to generate a >> video. Any direction? >> >> Angwenyi >> >> On 22 Dec 2015, at 12:10, Colin Cotter <colin.cotter(a)imperial.ac.uk> >> wrote: >> >> del output >> >> >> > > > -- > 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 > > >

2 1

Re: [firedrake] Animation
by Colin Cotter 22 Dec '15

22 Dec '15

Dear Angwenyi, The output should be one .pvd file, but many .vtu files, one per dump. Then, you open the .pvd file in Paraview, set up your visualisation and press the play button. all the best --cjc On 22 December 2015 at 12:06, Angwenyi David <kipkoej(a)gmail.com> wrote: > Dear Cotter, > > I did so, but got a single file. I thought I could get as many files as > timesteps; which files I could open all at once in Paraview to generate a > video. Any direction? > > Angwenyi > > On 22 Dec 2015, at 12:10, Colin Cotter <colin.cotter(a)imperial.ac.uk> > wrote: > > del output > > > -- http://www.imperial.ac.uk/people/colin.cotter www.cambridge.org/9781107663916

2 1

Re: [firedrake] Animation
by Andrew McRae 22 Dec '15

22 Dec '15

The code should look like this: output = File("waveEqn.pvd") output << phi while t < T: ... t += dt step += 1 if step % 10 == 0: output << phi On 22 December 2015 at 12:06, Angwenyi David <kipkoej(a)gmail.com> wrote: > Dear Cotter, > > I did so, but got a single file. I thought I could get as many files as > timesteps; which files I could open all at once in Paraview to generate a > video. Any direction? > > Angwenyi > > On 22 Dec 2015, at 12:10, Colin Cotter <colin.cotter(a)imperial.ac.uk> > wrote: > > del output > > >

1 0

Re: [firedrake] Animation
by Colin Cotter 22 Dec '15

22 Dec '15

Dear Angwenyi, You are overwriting the file in each step. To avoid this, you need to take the lines: ... File("waveEqn.pvd")<< phi ... output = File("waveEqn.pvd") outside of the time loop. all the best --cjc On 22 December 2015 at 11:05, Angwenyi David <kipkoej(a)gmail.com> wrote: > Dear all, > > I would like to get the animation on the linear wave equation in the > tutorial as a visualisation in Paraview. I used the following code to store > result files: > > ... phi -= dt / 2 * p > ... t += dt > ... if step % 10 == 0: > ... File("waveEqn.pvd")<< phi > ... output = File("waveEqn.pvd") > ... output << phi > ... del output > > I expected a number of files that would form snapshots for the animation, > but I got a single file. Could you please show me how to go about it. > > Regards. > > David > -- http://www.imperial.ac.uk/people/colin.cotter www.cambridge.org/9781107663916

2 1

Animation
by Angwenyi David 22 Dec '15

22 Dec '15

Dear all, I would like to get the animation on the linear wave equation in the tutorial as a visualisation in Paraview. I used the following code to store result files: ... phi -= dt / 2 * p ... t += dt ... if step % 10 == 0: ... File("waveEqn.pvd")<< phi ... output = File("waveEqn.pvd") ... output << phi ... del output I expected a number of files that would form snapshots for the animation, but I got a single file. Could you please show me how to go about it. Regards. David

1 0