June 2016 - firedrake - imperial.ac.uk

Saving txt files while running in parallel
by Floriane Gidel [RPG] 08 Jun '16

08 Jun '16

Dear all, I am running a Firedrake code in parallel with 4 cores, in which I save the maximal value of the amplitude of the wave at each time, with the command: max_eta = max(eta.dat.data) Eta_file.write('%-10s %-10s %-10s\n' % (t,'', eta_max)) But when opening the .txt file, I notice two issues: - the data are not saved continuously: time goes from 0 to 100 (let's say) and then starts from 90 again; - the values of max_eta are not the maximal values of eta, except at the beginning. So it looks like the function max takes the maximum value of eta only in one subdomain, so that after the wave has crossed the subdomain, the maximal value of eta goes back to the depth of rest value. How can I force this command to be applied on the full domain, even if the code is run in parallel? Thanks, Floriane

2 1

Gusto site up
by David Ham 07 Jun '16

07 Jun '16

Hi all, as requested, the Gusto site is up. If Colin provides me with a public key, I'll set things up so he can push updates Cheers, David

1 0

Fwd: [AMCG-Seminars] AMCG seminar on THURSDAY 9th June 4pm G38!
by Lawrence Mitchell 06 Jun '16

06 Jun '16

Dear all, this thursday Mark Adams will be giving an interesting talk on a new benchmarking metric for HPC systems, HPGMG. I'm sure it will be worthwhile attending! Lawrence > Begin forwarded message: > > From: Matthew Piggott <m.d.piggott(a)imperial.ac.uk> > Subject: Fwd: [AMCG-Seminars] AMCG seminar on THURSDAY 9th June 4pm G38! > Date: 6 June 2016 at 13:47:04 BST > To: Paul Kelly <p.kelly(a)imperial.ac.uk>, David Ham <david.ham(a)imperial.ac.uk>, "Mitchell, Lawrence" <lawrence.mitchell(a)imperial.ac.uk> > > You may want to advertise this to the appropriate lists. > Venue is G38 in the RSM building. > Best wishes > Matt > > > ---------- Forwarded message ---------- > From: Derrick, James <james.derrick10(a)imperial.ac.uk <mailto:james.derrick10@imperial.ac.uk>> > Date: 6 June 2016 at 11:46 > Subject: [AMCG-Seminars] AMCG seminar on THURSDAY 9th June 4pm G38! > To: amcg-seminars <amcg-seminars(a)imperial.ac.uk <mailto:amcg-seminars@imperial.ac.uk>> > > > Hi all, > > > This week's seminar will be on Thursday the 9th June (a day early this week) in G38, at the normal time of 4pm where Mark Adams will be giving a talk titled > > > "High Performance Geometric Multigrid: A Supercomputer Benchmark and Metric". > > > The abstract is below. > > Mark Adams and Sam Williams > > The HPGMG project (hpgmg.org <http://www.hpgmg.org/>) develops the extreme-scale benchmark and metric - HPGMG - which is a compact, stand alone, conceptually simple, benchmark, modeled in many ways on HPL. HPGMG solves a system of linear algebraic equations; the metric is equations solved per second (HPL maps equations to floating operations to generate the familiar flops/sec metric). HPGMG solves a finite volume discretization of the Laplacian and uses an appropriate solution algorithm - the non-iterative form of multigrid (full multigrid). HPGMG is designed for use as both a ranking metric and an aid to computer designers for developing architectures, and to centers for configuring systems that are well balanced with respect to general sparse applications. We publish the list twice each year (during ISC and SC). We present a preview of the ISC16 list, which shows some rank reordering relative to the Top500 list, and includes the first GPU (MPI+OpenMP+CUDA) results running on Titan and Piz Daint. > > > See you there! > > James & Roan > > > > > _______________________________________________ > AMCG and Friends Seminar Series Mailing List > AMCG-Seminars(a)imperial.ac.uk > Subscription Options: https://mailman.ic.ac.uk/mailman/listinfo/amcg-seminars > Webpage: http://amcg.ese.ic.ac.uk/index.php?title=Seminar_Series

1 0

PETSc Segmentation Violation
by Anna Kalogirou 06 Jun '16

06 Jun '16

Dear all, I am trying to run a code which returns the following error, before even entering the time loop: [0]PETSC ERROR: ------------------------------------------------------------------------ [0]PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation, probably memory access out of range [0]PETSC ERROR: Try option -start_in_debugger or -on_error_attach_debugger [0]PETSC ERROR: or see http://www.mcs.anl.gov/petsc/documentation/faq.html#valgrind [0]PETSC ERROR: or try http://valgrind.org on GNU/linux and Apple Mac OS X to find memory corruption errors [0]PETSC ERROR: configure using --with-debugging=yes, recompile, link, and run [0]PETSC ERROR: to get more information on the crash. -------------------------------------------------------------------------- MPI_ABORT was invoked on rank 0 in communicator MPI_COMM_WORLD with errorcode 59. NOTE: invoking MPI_ABORT causes Open MPI to kill all MPI processes. You may or may not see output from other processes, depending on exactly when Open MPI kills them. -------------------------------------------------------------------------- I tracked down the problem and I found that it is caused by the following: Matrix = assemble(u*v*dx).M.handle where u is a trial and v is a test function. Can any of you help with what might be the issue here? Thank you, Anna. -- Dr Anna Kalogirou Research Fellow School of Mathematics University of Leeds http://www1.maths.leeds.ac.uk/~matak/

5 9

Re: [firedrake] direct solver for mixed space runs out of memory
by Andrew McRae 03 Jun '16

03 Jun '16

Your FunctionSpace W has roughly 388*584*2 = 450k degrees of freedom. Your VectorFunctionSpace V has roughly 2*450k = 900k degrees of freedom. Then you construct a MixedFunctionSpace Q = W*W*V, which will have roughly 1.8 million degrees of freedom. LU causes fill-in, so (unless there's something very special about your PDE/discretisation) each row of the matrix will have O(sqrt(N)) [I think] non-zero entries which must be stored. And 1.8 million * O(sqrt(1.8 million)) * (8 bytes/double + 4 bytes/index) is probably more than 8 GB. Using MUMPS rather than the inbuilt LU might help? I.e., add the option "pc_factor_mat_solver_package": "mumps" to your solver_parameters dict. On 3 June 2016 at 23:15, Tianjiao Sun <tianjiao.sun14(a)imperial.ac.uk> wrote: > 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 >

2 1

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