- firedrake - imperial.ac.uk

direct solver for mixed space runs out of memory
by Tianjiao Sun 03 Jun '16

03 Jun '16

Hi team, I have the attached python code which is a mixed formulation on extruded mesh. It is possible that there are other bugs in my formulation, but I'm running into problems using the non-linear solver. My intent is to use direct solver in the non-linear iteration, and my parameters are: solver_parameters={'ksp_type': 'preonly', 'pc_type': 'lu', 'snes_atol':1e-10, 'ksp_atol':1e-10}) Then PETSc complains that I couldn't use direct solver for nested matrix. Google guides me to a post by Christian and Lawrence https://github.com/firedrakeproject/firedrake/issues/431 So I added nest=False in solve(), but then the system run out of memory: [0] Out of memory. Allocated: 0, Used by process: 5061259264 [0] Memory requested 8589934188 I guess it's not quite likely to use 8 GB of RAM for this problem size (388 x 584 x 2 layers). Lawrence did say that "it doesn't work with vector function spaces" earlier in the post, I'm not sure if that is related? My mixed space is [W, W, V] where W = DG(0) x CG(1), V is vector space of dim=2. Many thanks, -TJ

1 0

Dimension of the gradient of a vector
by Floriane Gidel [RPG] 02 Jun '16

02 Jun '16

Dear all, I am trying to understand how the gradient operator works on Firedrake and there is something I misunderstand. As far as I understood, if I define (in the HORIZONTAL plane): V = FunctionSpace(mesh,"CG", 1) and W = VectorFunctionSpace(mesh,"CG", 1, dim=N) such that : psi_1 = Function (V) psi_N = Function(W) (vector of dim (1,N)) then, grad(psi_1) is a rank 1 tensor of 2 rows or columns (does not matter), grad(psi_N) is a rank 2 tensor of dim (N, 2) 1) Is that correct ? 2) Then, I understand that we can compute grad(psi_N)*grad(psi_1) but not grad(psi_1)*grad(psi_N), because the dimensions agree in the first case but not in the second one. However, when taking the transpose of grad(psi_N), I don't understand what happens. The product: grad(psi_1)*(grad(psi_N)).T cannot be computed, and the product (grad(psi_N)).T * grad(psi_1) can be computed, while I would expect the contrary (since I expect (grad(psi_N)).T to be of dimension (2,N) ). Moreover, if I print the expression of (grad(psi_N)).T*grad(psi_1), I get : { A | A_{i_8} = sum_{i_9} ((grad(w_6))^T)[i_8, i_9] * (grad(w_4))[i_9] } which is actually A | A_i = sum_{j=1}^{j=2} (grad(psi_N))^T(i,j) * d psi_1/dxj But if grad(psi_N) has dimension (N,2), then (grad(psi_N))^T has dimension (2,N) and therefore the expression above would not work if N>2. Can someone explain me how the gradient of psi_N and its transpose are evaluated, and what are their rank and dimension ? Is there a command to print these dimensions ? Thanks and best regards, Floriane

3 3

FEniCS implementation works but Firedrake does not, why?
by Justin Chang 01 Jun '16

01 Jun '16

Hi all, So we have been attempting to convert a FEniCS implementation of a semi-linear diffusion code into Firedrake. Mainly because we want to employ PETSc/TAO's optimization routines, which FEniCS does not let us do. However, the Firedrake implementation is not working, namely our consistent Newton-Raphson approach. Attached is the FEniCS code (P3_Galerkin_NR.py) and the firedrake code (NR_Nonlinear_poisson.py). For FEniCS you may just run the code as "python P3_Galerkin_NR.py" but for the firedrake code you must run as: python NR_Nonlinear_poisson.py 50 50 0 For the FEniCS code, this is my solver output: >> python P3_Galerkin_NR.py Calling FFC just-in-time (JIT) compiler, this may take some time. Calling FFC just-in-time (JIT) compiler, this may take some time. Solving linear variational problem. iter=1: norm=1 Solving linear variational problem. iter=2: norm=0.0731224 Solving linear variational problem. iter=3: norm=0.00217701 Solving linear variational problem. iter=4: norm=1.64398e-06 Solving linear variational problem. iter=5: norm=8.89289e-13 but for Firedrake, i have this: >> python NR_Nonlinear_poisson.py 50 50 0 COFFEE finished in 0.00142097 seconds (flops: 0 -> 0) COFFEE finished in 0.00103188 seconds (flops: 0 -> 0) COFFEE finished in 0.00173187 seconds (flops: 0 -> 0) COFFEE finished in 0.00173306 seconds (flops: 0 -> 0) COFFEE finished in 0.0015831 seconds (flops: 0 -> 0) COFFEE finished in 0.000972033 seconds (flops: 0 -> 0) COFFEE finished in 0.00118279 seconds (flops: 2 -> 2) COFFEE finished in 0.000946999 seconds (flops: 0 -> 0) COFFEE finished in 0.00165105 seconds (flops: 1 -> 1) Error norm: 2.193e-01 Error norm: 4.697e-02 Error norm: 4.830e-02 Error norm: 8.437e-02 Error norm: 2.740e-01 Error norm: 2.872e+00 Error norm: 8.993e+02 Error norm: 1.912e+10 Error norm: 1.919e+32 Error norm: 1.981e+98 Traceback (most recent call last): File "NR_Nonlinear_poisson.py", line 119, in <module> solver.solve(u_k1,b) File "/home/justin/Software/firedrake/src/firedrake/firedrake/linear_solver.py", line 153, in solve raise RuntimeError("LinearSolver failed to converge after %d iterations with reason: %s", self.ksp.getIterationNumber(), solving_utils.KSPReasons[r]) RuntimeError: ('LinearSolver failed to converge after %d iterations with reason: %s', 0, 'DIVERGED_NANORINF') I strongly believe there are no inconsistencies between our FEniCS and Firedrake codes (in terms of what we want to solve) but for some reason the latter blows up. However, if I run the optimization version of the Firedrake code, my solver converges quadratically as expected: >> python NR_Nonlinear_poisson.py 50 50 1 COFFEE finished in 0.00143099 seconds (flops: 0 -> 0) COFFEE finished in 0.00164199 seconds (flops: 0 -> 0) COFFEE finished in 0.00169611 seconds (flops: 0 -> 0) COFFEE finished in 0.00169706 seconds (flops: 0 -> 0) COFFEE finished in 0.00153899 seconds (flops: 0 -> 0) COFFEE finished in 0.000952005 seconds (flops: 0 -> 0) COFFEE finished in 0.000900984 seconds (flops: 1 -> 1) COFFEE finished in 0.000989914 seconds (flops: 0 -> 0) Error norm: 2.193e-01 Error norm: 1.626e-02 Error norm: 4.870e-04 Error norm: 9.332e-07 Error norm: 0.000e+00 When I compare the pvd plots between the FEniCS and Firedrake (with optimization), Firedrake seems correct qualitatively (i.e., the negative concentrations are gone). But I am confused why the non-optimization implementation of firedrake does not converge whereas the FEniCS one does. Any thoughts? Thanks! Justin

2 4

Re: [firedrake] FEniCS implementation works but Firedrake does not, why?
by Colin Cotter 01 Jun '16

01 Jun '16

Heh, ok. It's just that it quickly alerted me to some situations where COFFEE was actually producing code with more FLOPS! --cjc On 1 June 2016 at 18:22, Ham, David A <david.ham(a)imperial.ac.uk> wrote: > On the contrary. Please file a bug against Coffee. If enough people > complain Fabio might shut the bloody thing up! ;) > > David > On Wed, 1 Jun 2016 at 18:19, Cotter, Colin J <colin.cotter(a)imperial.ac.uk> > wrote: > >> Hi Justin, >> We'd prefer you not to disable these messages since they help to notify >> us if there are problems with the COFFEE optimisations. >> >> all the best >> --cjc >> >> On 1 June 2016 at 18:12, Justin Chang <jychang48(a)gmail.com> wrote: >> >>> 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 >>>>> >>>>> >>>> >>> >> >> >> -- >> 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] FEniCS implementation works but Firedrake does not, why?
by David Ham 01 Jun '16

01 Jun '16

On the contrary. Please file a bug against Coffee. If enough people complain Fabio might shut the bloody thing up! ;) David On Wed, 1 Jun 2016 at 18:19, Cotter, Colin J <colin.cotter(a)imperial.ac.uk> wrote: > Hi Justin, > We'd prefer you not to disable these messages since they help to notify > us if there are problems with the COFFEE optimisations. > > all the best > --cjc > > On 1 June 2016 at 18:12, Justin Chang <jychang48(a)gmail.com> wrote: > >> 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 >>>> >>>> >>> >> > > > -- > http://www.imperial.ac.uk/people/colin.cotter > > www.cambridge.org/9781107663916 > > >

1 0

Re: [firedrake] FEniCS implementation works but Firedrake does not, why?
by Colin Cotter 01 Jun '16

01 Jun '16

Hi Justin, We'd prefer you not to disable these messages since they help to notify us if there are problems with the COFFEE optimisations. all the best --cjc On 1 June 2016 at 18:12, Justin Chang <jychang48(a)gmail.com> wrote: > 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 >>> >>> >> > -- http://www.imperial.ac.uk/people/colin.cotter www.cambridge.org/9781107663916

2 1

Re: [firedrake] FEniCS implementation works but Firedrake does not, why?
by David Ham 01 Jun '16

01 Jun '16

I confirm I can reproduce the firedrake results. I'll try to have a poke to see if I can work out what is going wrong. David On Wed, 1 Jun 2016 at 09:52 Mitchell, Lawrence < lawrence.mitchell(a)imperial.ac.uk> wrote: > > > On 01/06/16 09:44, Fabio Luporini wrote: > > Hi Justin, > > > > since you're talking about possible issues with > > optimised/non-optimised Firedrake, I checked if this had something to > > do with COFFEE, but it really seems it doesn't, so I'm not sure what's > > going on. > > Different type of optimisation! > > Lawrence > >

1 0

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

01 Jun '16

Hi Justin, since you're talking about possible issues with optimised/non-optimised Firedrake, I checked if this had something to do with COFFEE, but it really seems it doesn't, so I'm not sure what's going on. Maybe the others can help. -- Fabio 2016-06-01 9:44 GMT+02:00 Justin Chang <jychang48(a)gmail.com>: > Hi all, > > So we have been attempting to convert a FEniCS implementation of a > semi-linear diffusion code into Firedrake. Mainly because we want to employ > PETSc/TAO's optimization routines, which FEniCS does not let us do. > However, the Firedrake implementation is not working, namely our consistent > Newton-Raphson approach. > > Attached is the FEniCS code (P3_Galerkin_NR.py) and the firedrake code > (NR_Nonlinear_poisson.py). For FEniCS you may just run the code as "python > P3_Galerkin_NR.py" but for the firedrake code you must run as: > > python NR_Nonlinear_poisson.py 50 50 0 > > For the FEniCS code, this is my solver output: > > >> python P3_Galerkin_NR.py > Calling FFC just-in-time (JIT) compiler, this may take some time. > Calling FFC just-in-time (JIT) compiler, this may take some time. > Solving linear variational problem. > iter=1: norm=1 > Solving linear variational problem. > iter=2: norm=0.0731224 > Solving linear variational problem. > iter=3: norm=0.00217701 > Solving linear variational problem. > iter=4: norm=1.64398e-06 > Solving linear variational problem. > iter=5: norm=8.89289e-13 > > but for Firedrake, i have this: > > >> python NR_Nonlinear_poisson.py 50 50 0 > COFFEE finished in 0.00142097 seconds (flops: 0 -> 0) > COFFEE finished in 0.00103188 seconds (flops: 0 -> 0) > COFFEE finished in 0.00173187 seconds (flops: 0 -> 0) > COFFEE finished in 0.00173306 seconds (flops: 0 -> 0) > COFFEE finished in 0.0015831 seconds (flops: 0 -> 0) > COFFEE finished in 0.000972033 seconds (flops: 0 -> 0) > COFFEE finished in 0.00118279 seconds (flops: 2 -> 2) > COFFEE finished in 0.000946999 seconds (flops: 0 -> 0) > COFFEE finished in 0.00165105 seconds (flops: 1 -> 1) > Error norm: 2.193e-01 > Error norm: 4.697e-02 > Error norm: 4.830e-02 > Error norm: 8.437e-02 > Error norm: 2.740e-01 > Error norm: 2.872e+00 > Error norm: 8.993e+02 > Error norm: 1.912e+10 > Error norm: 1.919e+32 > Error norm: 1.981e+98 > Traceback (most recent call last): > File "NR_Nonlinear_poisson.py", line 119, in <module> > solver.solve(u_k1,b) > File > "/home/justin/Software/firedrake/src/firedrake/firedrake/linear_solver.py", > line 153, in solve > raise RuntimeError("LinearSolver failed to converge after %d > iterations with reason: %s", self.ksp.getIterationNumber(), > solving_utils.KSPReasons[r]) > RuntimeError: ('LinearSolver failed to converge after %d iterations with > reason: %s', 0, 'DIVERGED_NANORINF') > > > I strongly believe there are no inconsistencies between our FEniCS and > Firedrake codes (in terms of what we want to solve) but for some reason the > latter blows up. However, if I run the optimization version of the > Firedrake code, my solver converges quadratically as expected: > > >> python NR_Nonlinear_poisson.py 50 50 1 > COFFEE finished in 0.00143099 seconds (flops: 0 -> 0) > COFFEE finished in 0.00164199 seconds (flops: 0 -> 0) > COFFEE finished in 0.00169611 seconds (flops: 0 -> 0) > COFFEE finished in 0.00169706 seconds (flops: 0 -> 0) > COFFEE finished in 0.00153899 seconds (flops: 0 -> 0) > COFFEE finished in 0.000952005 seconds (flops: 0 -> 0) > COFFEE finished in 0.000900984 seconds (flops: 1 -> 1) > COFFEE finished in 0.000989914 seconds (flops: 0 -> 0) > Error norm: 2.193e-01 > Error norm: 1.626e-02 > Error norm: 4.870e-04 > Error norm: 9.332e-07 > Error norm: 0.000e+00 > > When I compare the pvd plots between the FEniCS and Firedrake (with > optimization), Firedrake seems correct qualitatively (i.e., the negative > concentrations are gone). But I am confused why the non-optimization > implementation of firedrake does not converge whereas the FEniCS one does. > > Any thoughts? > > Thanks! > Justin >

2 1

Detecting NaNs
by Gregory, Alastair C A 31 May '16

31 May '16

Hi Firedrakers, I have two functions, v(x,y) and h(x,y). I find a quantity u(x,y)=v(x,y)/h(x,y) to solve a PDE. Now, I'm struggling with using u(x,y) when h(x,y) is equal to 0 at some cells. I get DIVERGED_NAN errors as I'm trying to divide by 0. Now what I have done is to build a kernel that allows me to set the cells that will be dividing by 0, to 0, and otherwise proceed as normal on the cells that won't be dividing by 0. This arises since v is momentum and h is depth, and thus u is velocity. With 0 depth, I would like 0 velocity. However, sometimes I divide by a restriction to a function and thus can't use a kernel. I was wondering if there was a clean way (potentially UFL way) to set cells of functions from NANS to 0, or anything similar that would help me? Many Thanks, Alastair Gregory ------------------------------------------------------------------------------------------ Research Assistant (Maths Foresees Feasibility Study) Imperial College London Room 759 Huxley Building, South Kensington (Tel: 07794 243913) | (Email: a.gregory14(a)imperial.ac.uk) -------------------------------------------------------------------------------------------

2 1

Two coupled problems
by Daniel Ruprecht 31 May '16

31 May '16

Dear Firedrakers, I am trying to solve two coupled problems on different domains (essentially a small box next to a larger box), where the solution at the boundary from one problem determines the Neumann boundary condition for the second problem at the interface. I am not sure how to efficiently use the solution of the first problem to prescribe the boundary condition for the second. My idea was to derive from the Expression class and overwrite the eval function including a couple of if statements to make sure the coordinates are in the region where both meshes actually match. Is this the correct way to do it (it seems quite slow, probably because the derived expression can't be compiled in C) or is there a cleverer way to do this? Cheers, Daniel -- Dr. Daniel Ruprecht University Academic Fellow School of Mechanical Engineering, Room 4.50 University of Leeds Leeds LS2 9JT, UK Email: d.ruprecht(a)leeds.ac.uk Phone: +44 (0)113 343 22 01

2 6