- firedrake - imperial.ac.uk

optimisation of the running time
by Floriane Gidel [RPG] 12 Jul '16

12 Jul '16

Dear all, I would like some advice to decrease the running time in my code 3D_NL.py . Details are given in the attached pdf and below is a quick summary: In the attached files, I run a 3D nonlinear code for potential flow water wave equations. The vertical and horizontal space discretisation are separated, so that Firedrake solves the equations in the horizontal plane only, and the vertical dependency is involved through matrices computed by evaluating the z-dependent integrals with the functions defined in the file 'discr_sigma.py'. Basically, the potential flow phi(x,y,z,t) is defined as phi(x,y,t)phi(z), where phi(z) is discretised as a n^th order Lagrange expansion in one element: phi(z_i) = product_{k=0}^n (z-z_k)/(z_i-z_k) for k!= i There are therefore 3 unknowns, which are : - h (the total depth of water) defined in the horizontal plane - phi_s (the potential flow at the surface) defined in the horizontal plane - hat_psi a vector function containing all the phi_i for i = 1:n , which are the values of phi in the vertical nodes (excluding the surface ones), where n is the order of the Lagrange expansion. Each phi_i is defined in the horizontal plane. The time scheme is defined like this : step 1 : update phi_s and hat_psi simultaneously (that is, sum of n+1 weak formulations) step 2 : update h and hat_psi simultaneously (that is, sum of n+1 weak formulations) step 3 : update phi_s (one weak formulation). (those steps are written in the pdf file as well) At the moment, the equations are solved in a small domain with resolution 0.1, but it's running very slowly. Running in parallel does not improve it since the number of nodes is too small. Is there another way to increase the simulation speed? Maybe by changing the iteration criteria in the solver ? I hope this is clear enough, otherwise I can send you more details. Best regards, Floriane

3 5

Re: [firedrake] Using Firedrake to model optical fibres
by David Ham 05 Jul '16

05 Jul '16

Hi Ian, Thanks for your interest in Firedrake. You could definitely solve Maxwell's equations in real space in Firedrake, even using the right exterior calculus discrete function spaces. It will also work in parallel for free. I'm less sure whether we can do Fourier space. I think I'd need to see the formulation. In particular, we can't currently solve PDEs over the complex numbers - although there are some work arounds for this. (Complex is definitely on the wish list but it will only happen if I manage to attract the right masters student to work on it). Regards, David On Tue, 5 Jul 2016 at 14:35 Ian Thompson <I.R.Thompson(a)bath.ac.uk> wrote: > Hello all, > > I saw Firedrake at the Bath HPC day, I work in physics and we do a lot of > fibre optic research. This includes modelling Maxwell’s equations across a > 2D mesh that describes a cross-section of a fibre. Traditionally we solve > this using a finite element method and solving in Fourier space. Is > Firedrake suitable for solving such problems? > Looking at the documentation Firedrake can load meshes from other > programs, to solve in Fourier space would it be sufficient to define the > equations in Fourier space and load a mesh that corresponds to the Fourier > transform of our real space structure. Or can Firedrake support the > transform itself? Would Firedrake be capable of solving the equations in > real-space? In principle there is no problem with solving in real-space, > its just been neater in the past. > Currently we use COMSOL which is expensive and slow due to a lot of > unnecessary (for our purposes at least) overhead and lack of parallel > licensing. Sorry for so many questions but I am quietly excited at using a > neater, nicer solver. > > Cheers > Ian > _______________________________________________ > firedrake mailing list > firedrake(a)imperial.ac.uk > https://mailman.ic.ac.uk/mailman/listinfo/firedrake >

1 0

Using Firedrake to model optical fibres
by Ian Thompson 05 Jul '16

05 Jul '16

Hello all, I saw Firedrake at the Bath HPC day, I work in physics and we do a lot of fibre optic research. This includes modelling Maxwell’s equations across a 2D mesh that describes a cross-section of a fibre. Traditionally we solve this using a finite element method and solving in Fourier space. Is Firedrake suitable for solving such problems? Looking at the documentation Firedrake can load meshes from other programs, to solve in Fourier space would it be sufficient to define the equations in Fourier space and load a mesh that corresponds to the Fourier transform of our real space structure. Or can Firedrake support the transform itself? Would Firedrake be capable of solving the equations in real-space? In principle there is no problem with solving in real-space, its just been neater in the past. Currently we use COMSOL which is expensive and slow due to a lot of unnecessary (for our purposes at least) overhead and lack of parallel licensing. Sorry for so many questions but I am quietly excited at using a neater, nicer solver. Cheers Ian

1 0

Re: [firedrake] making a broken element from a Function Space
by Colin Cotter 05 Jul '16

05 Jul '16

OK, this seems to work. I made a PR in FIAT. cheers --cjc On 5 July 2016 at 11:01, Cotter, Colin J <colin.cotter(a)imperial.ac.uk> wrote: > Ah, I assumed it was being tested in the hybridisation test but then I > remembered that one doesn't do extruded elements. > > cheers > --cjc > > On 5 July 2016 at 10:27, Mitchell, Lawrence < > lawrence.mitchell(a)imperial.ac.uk> wrote: > >> >> > On 5 Jul 2016, at 10:20, Colin Cotter <colin.cotter(a)imperial.ac.uk> >> wrote: >> > >> > Yes, it then just hit another use of ref_el lower down. I can fix this >> of course, but I wonder if I am doing something weird since this is just >> the usual usage of BrokenElement and hasn't been an issue in other examples. >> >> >> They're not used in any tests on extruded elements. Which is probably >> why this issue hasn't appeared. >> >> Lawrence >> > > > > -- > http://www.imperial.ac.uk/people/colin.cotter > > www.cambridge.org/9781107663916 > > > -- http://www.imperial.ac.uk/people/colin.cotter www.cambridge.org/9781107663916

1 0

Re: [firedrake] making a broken element from a Function Space
by Colin Cotter 05 Jul '16

05 Jul '16

Ah, I assumed it was being tested in the hybridisation test but then I remembered that one doesn't do extruded elements. cheers --cjc On 5 July 2016 at 10:27, Mitchell, Lawrence < lawrence.mitchell(a)imperial.ac.uk> wrote: > > > On 5 Jul 2016, at 10:20, Colin Cotter <colin.cotter(a)imperial.ac.uk> > wrote: > > > > Yes, it then just hit another use of ref_el lower down. I can fix this > of course, but I wonder if I am doing something weird since this is just > the usual usage of BrokenElement and hasn't been an issue in other examples. > > > They're not used in any tests on extruded elements. Which is probably why > this issue hasn't appeared. > > Lawrence > -- http://www.imperial.ac.uk/people/colin.cotter www.cambridge.org/9781107663916

1 0

Re: [firedrake] making a broken element from a Function Space
by Colin Cotter 05 Jul '16

05 Jul '16

Yes, it then just hit another use of ref_el lower down. I can fix this of course, but I wonder if I am doing something weird since this is just the usual usage of BrokenElement and hasn't been an issue in other examples. --cjc On 5 July 2016 at 10:00, Mitchell, Lawrence < lawrence.mitchell(a)imperial.ac.uk> wrote: > > > On 5 Jul 2016, at 09:41, Colin Cotter <colin.cotter(a)imperial.ac.uk> > wrote: > > > > Actually, I'm now less sure where the problem is. > > > > If I do this, in the state object > > ... > > self.V_elt[2] = V2t_elt > > ... > > > > and then later, > > Vtdg_elt = BrokenElement(state.V_elt[2]) > > Vtdg = FunctionSpace(mesh, Vtdg_elt) > > > > then I get the same error. > > Did you try my suggestion? > > Lawrence > -- http://www.imperial.ac.uk/people/colin.cotter www.cambridge.org/9781107663916

2 1

Re: [firedrake] making a broken element from a Function Space
by Colin Cotter 05 Jul '16

05 Jul '16

Actually, I'm now less sure where the problem is. If I do this, in the state object ... self.V_elt[2] = V2t_elt ... and then later, Vtdg_elt = BrokenElement(state.V_elt[2]) Vtdg = FunctionSpace(mesh, Vtdg_elt) then I get the same error. --cjc On 5 July 2016 at 09:31, Mitchell, Lawrence < lawrence.mitchell(a)imperial.ac.uk> wrote: > > > On 5 Jul 2016, at 09:23, Colin Cotter <colin.cotter(a)imperial.ac.uk> > wrote: > > > > Dear Firedrakers, > > I have a bit of code where I want to construct a broken function space > corresponding to a given function space. > > > > My code is: > > Vdg_elt = BrokenElement(V.ufl_element()) > > Vdg = FunctionSpace(state.mesh, Vdg_elt) > > > > This returns an error: > > File "/home/cjc1/repositories/dcore/gusto/advection.py", line 199, in > __init__ > > Vdg = FunctionSpace(state.mesh, Vdg_elt) > > File "<decorator-gen-291>", line 2, in FunctionSpace > > File "/home/cjc1/firedrake/src/PyOP2/pyop2/profiling.py", line 59, in > wrapper > > return f(*args, **kwargs) > > File "/home/cjc1/firedrake/src/firedrake/firedrake/functionspace.py", > line 131, in FunctionSpace > > new = impl.FunctionSpace(topology, element, name=name) > > File > "/home/cjc1/firedrake/src/firedrake/firedrake/functionspaceimpl.py", line > 245, in __init__ > > sdata = get_shared_data(mesh, fiat_element) > > File > "/home/cjc1/firedrake/src/firedrake/firedrake/functionspacedata.py", line > 558, in get_shared_data > > return FunctionSpaceData(mesh, fiat_element) > > File > "/home/cjc1/firedrake/src/firedrake/firedrake/functionspacedata.py", line > 310, in __init__ > > self.bt_masks = get_bt_masks(mesh, edofs_key, fiat_element) > > File "<decorator-gen-281>", line 2, in get_bt_masks > > File > "/home/cjc1/firedrake/src/firedrake/firedrake/functionspacedata.py", line > 56, in cached > > result = f(mesh, key, *args, **kwargs) > > File > "/home/cjc1/firedrake/src/firedrake/firedrake/functionspacedata.py", line > 205, in get_bt_masks > > facet_dofs = horiz_facet_support_dofs(fiat_element) > > File "/home/cjc1/firedrake/src/fiat/FIAT/tensor_product.py", line 355, > in horiz_facet_support_dofs > > assert isinstance(elem.ref_el, TensorProductCell) > > AttributeError: DiscontinuousElement instance has no attribute 'ref_el' > > > > What is the right way to do this? > > > Patch FIAT/tensor_product.py to use elem.get_reference_element() rather > than using internal implementation details. > > Lawrence > -- http://www.imperial.ac.uk/people/colin.cotter www.cambridge.org/9781107663916

2 1

making a broken element from a Function Space
by Colin Cotter 05 Jul '16

05 Jul '16

Dear Firedrakers, I have a bit of code where I want to construct a broken function space corresponding to a given function space. My code is: Vdg_elt = BrokenElement(V.ufl_element()) Vdg = FunctionSpace(state.mesh, Vdg_elt) This returns an error: File "/home/cjc1/repositories/dcore/gusto/advection.py", line 199, in __init__ Vdg = FunctionSpace(state.mesh, Vdg_elt) File "<decorator-gen-291>", line 2, in FunctionSpace File "/home/cjc1/firedrake/src/PyOP2/pyop2/profiling.py", line 59, in wrapper return f(*args, **kwargs) File "/home/cjc1/firedrake/src/firedrake/firedrake/functionspace.py", line 131, in FunctionSpace new = impl.FunctionSpace(topology, element, name=name) File "/home/cjc1/firedrake/src/firedrake/firedrake/functionspaceimpl.py", line 245, in __init__ sdata = get_shared_data(mesh, fiat_element) File "/home/cjc1/firedrake/src/firedrake/firedrake/functionspacedata.py", line 558, in get_shared_data return FunctionSpaceData(mesh, fiat_element) File "/home/cjc1/firedrake/src/firedrake/firedrake/functionspacedata.py", line 310, in __init__ self.bt_masks = get_bt_masks(mesh, edofs_key, fiat_element) File "<decorator-gen-281>", line 2, in get_bt_masks File "/home/cjc1/firedrake/src/firedrake/firedrake/functionspacedata.py", line 56, in cached result = f(mesh, key, *args, **kwargs) File "/home/cjc1/firedrake/src/firedrake/firedrake/functionspacedata.py", line 205, in get_bt_masks facet_dofs = horiz_facet_support_dofs(fiat_element) File "/home/cjc1/firedrake/src/fiat/FIAT/tensor_product.py", line 355, in horiz_facet_support_dofs assert isinstance(elem.ref_el, TensorProductCell) AttributeError: DiscontinuousElement instance has no attribute 'ref_el' What is the right way to do this? all the best --cjc -- http://www.imperial.ac.uk/people/colin.cotter www.cambridge.org/9781107663916

2 1

Assembling Jacobian within TAO
by Justin Chang 30 Jun '16

30 Jun '16

Hi all, So I am attempting to use TAO's mixed complementarity solvers and it basically requires a formGradient/Constraints routine and a formJacobian routine. I am solving a variational inequality via TAO and want the linear operator A to be my jacobian. However, I don't think it was correctly done because my semi smooth linear solver does not converge (but it should!). Attached is the code. Line 108 is what TaoSetJacobianRoutine() calls. Basically I want to assemble the linear operator into the matrix petsc_J. Is what I have the correct implementation? run the code as: > python MCP_poisson.py 100 100 2 Thanks, Justin

2 2

Re: [firedrake] FEniCS implementation works but Firedrake does not, why?
by Fabio Luporini 29 Jun '16

29 Jun '16

Hi Justin, I should be able to come up with a non-verbose mode in the next few days, thanks -- Fabio 2016-06-01 19:12 GMT+02:00 Justin Chang <jychang48(a)gmail.com>: > A couple more related questions: > > 1) Every time I run this (or any firedrake code) I get these notifications: > > COFFEE finished in XXX seconds (flops: X -> X) > > Is there a way to disable these? > > 2) Also, if I use quadrilateral elements, i also get this notification (I > believe it happens during the errornorm() call) > > Discontinuous Lagrange element requested on quadrilateral, creating DQ > element. > > Is there a way to also disable this? > > Thanks, > Justin > > On Wed, Jun 1, 2016 at 11:56 AM, Justin Chang <jychang48(a)gmail.com> wrote: > >> Ah this works wonderfully, thanks Lawrence! >> >> On Wed, Jun 1, 2016 at 4:23 AM, Lawrence Mitchell < >> lawrence.mitchell(a)imperial.ac.uk> wrote: >> >>> Hi Justin, >>> >>> On 01/06/16 08:44, Justin Chang wrote: >>> >>> ... >>> >>> > #===================== >>> > # Medium properties >>> > #===================== >>> > D = Function(Q) >>> > alpha = Constant(0.0) >>> > d1 = 1.0 >>> > d2 = 0.0001 >>> > theta = pi/6.0 >>> > co = cos(theta) >>> > si = sin(theta) >>> > >>> D.interpolate(Expression(("co*co*d1+si*si*d2","-co*si*(d1-d2)","-co*si*(d1-d2)","si*si*d1+co*co*d2"),co=co,si=si,d1=d1,d2=d2)) >>> >>> >>> FWIW, I recommend moving to using interpolation of UFL expressions for >>> these kind of initialisations (you get, amongst other things) better >>> error messages when mistyping something. See >>> http://firedrakeproject.org/interpolation.html#ufl-expressions >>> >>> > #===================== >>> > # Volumetric source >>> > #===================== >>> > def f(u): >>> > return Constant(10)*u*u*u >>> > >>> > def df(u): >>> > return Constant(30)*u*u >>> > >>> > #===================== >>> > # Variational form >>> > #===================== >>> > a = alpha * inner(u,v) * dx + inner(D*nabla_grad(u), nabla_grad(v))*dx >>> + inner(df(u_k) * u,v)*dx >>> >>> So a is a bilinear form, but it depends on the value of the previous >>> solution. >>> >>> > A = assemble(a,bcs=bcs) >>> > solver = LinearSolver(A,options_prefix="solver_") >>> >>> Here you assemble A and build a LinearSolver. >>> >>> > #======================== >>> > # Optimization >>> > #======================== >>> > lb = Function(V) >>> > ub = Function(V) >>> > ub.assign(1000) >>> > taoSolver = PETSc.TAO().create(PETSc.COMM_WORLD) >>> > with lb.dat.vec_ro as lb_vec, ub.dat.vec as ub_vec: >>> > taoSolver.setVariableBounds(lb_vec,ub_vec) >>> > >>> > def ObjGrad(tao, petsc_x, petsc_g): >>> > with b.dat.vec as b_vec: >>> > A.M.handle.mult(petsc_x,petsc_g) >>> > xtHx = petsc_x.dot(petsc_g) >>> > xtf = petsc_x.dot(b_vec) >>> > petsc_g.axpy(-1.0,b_vec) >>> > return 0.5*xtHx - xtf >>> > >>> > taoSolver.setObjectiveGradient(ObjGrad) >>> > taoSolver.setType(PETSc.TAO.Type.BLMVM) >>> > >>> > #======================== >>> > # Picard iteration >>> > #======================== >>> > while it < maxit and eps > tol: >>> > A = assemble(a,bcs=bcs) >>> >>> Here you reassemble a and assign it to a new variable "A". >>> >>> > # Standard solver >>> > if opt == 0: >>> > b = assemble(L,bcs=bcs) >>> > solver.solve(u_k1,b) >>> >>> This is never seen inside the solver, so this solver always runs with >>> the originally assembled A (with the initial value for u_k). >>> >>> > # Optimization solver >>> > else: >>> > b = assemble(L) >>> > b -= rhs_opt >>> > for bc in bcs: >>> > bc.apply(b) >>> > with u_k1.dat.vec as u_k1_vec: >>> > taoSolver.solve(u_k1_vec) >>> g >>> I'm actually slightly surprised that this one works, but that appears >>> to be because I don't understand how python's scoping rules work. >>> >>> The LinearSolver object isn't designed to work if your jacobian >>> changes hence the problem. >>> >>> I recommend you use a LinearVariationalSolver, but you need to be >>> careful, since by default it doesn't rebuild the jacobian every solve. >>> >>> If I do: >>> >>> problem = LinearVariationalProblem(a, L, u_k1, bcs=bcs, >>> constant_jacobian=False) >>> solver = LinearVariationalSolver(problem, options_prefix="solver_") >>> >>> if opt == 0: >>> solver.solve() >>> >>> >>> Then I get: >>> >>> python foo.py 50 50 0 >>> >>> Error norm: 2.193e-01 >>> Error norm: 1.705e-02 >>> Error norm: 5.105e-04 >>> Error norm: 4.980e-07 >>> Error norm: 5.332e-13 >>> >>> We could add the option to rebuild the jacobian inside a LinearSolver >>> (maybe if you've called something like solver.reset_operator() first), >>> but it's not available right now. >>> >>> I've attached a modified version that gets quadratic convergence in >>> both cases. I've also modified the taoSolver setup slightly so that >>> you reuse the same storage for the operator and RHS (rather than >>> building a new one every time). >>> >>> Cheers, >>> >>> Lawrence >>> >>> _______________________________________________ >>> firedrake mailing list >>> firedrake(a)imperial.ac.uk >>> https://mailman.ic.ac.uk/mailman/listinfo/firedrake >>> >>> >> >

3 4