Lawrence, Thanks for all of your help. Here’s what I ended up with: —————— def ObjGrad_D(tao, petsc_x, petsc_g): with b_D_nobc.dat.vec as b_D_vec: A_D_nobc.M.handle.mult(petsc_x,petsc_g) xtHx = petsc_x.dot(petsc_g) xtf = petsc_x.dot(b_D_vec) petsc_g.axpy(-1.0,b_D_vec) arr = petsc_g.array for bc in bc_D: nodes = bc.nodes[bc.nodes < bc.function_space().dof_dset.size] if isinstance(bc.function_arg, Function): with bc.function_arg.dat.vec_ro as bc_vec: arr[nodes] = bc_vec.array_r[nodes] else: arr[nodes] = float(bc.function_arg) return 0.5*xtHx - xtf ——————— However it seems something is still off about the boundary conditions I am applying. I have attached three different scenarios for my two solutions psi_A and psi_B. The boundary conditions are: psi_A: set to 1 at x == 0.0 and 0.0 <= y <= 0.5, 0 everywhere else psi_B: set to 1 at x == 0.0 and 0.5 < y <= 1.0, 0 everywhere else - no_opt_psiA/B.png illustrates solving the standard diffusion without any optimization. As far as I know, these solutions are “correct” in terms of the boundary conditions. - opt_wo_bounds_psiA/B.png illustrates solving the standard diffusion with TAO’s optimization but without any variable bounds. I would expect this to have the same solution as the no_opt cases but it seems to me the boundary conditions are messed up; face ID’s 3 and 4 (i.e., top and bottom) should have homogeneous bcs but it seems they are not registering. - opt_w_bounds_psiA/B.png now includes the variable bounds set by these lines: lb = Function(W) ub = Function(W) lb.assign(0) ub.assign(1) tao = PETSc.TAO().create(MPI.COMM_WORLD) with lb.dat.vec as lb_vec, ub.dat.vec as ub_vec: tao.setVariableBounds(lb_vec,ub_vec) The boundary conditions are also screwed up at the top and bottom faces. In this case, psi_B (corresponding to W.sub(1)) is even more screwed up - the non-homogeneous part of psi_B isn’t even registering (hence concentrations of < 1). Replacing the RHS of the arr[nodes] assignments with 0 didn’t change the solution either. Any idea what’s going on here? Thanks, Justin

On Sep 14, 2015, at 5:02 AM, Lawrence Mitchell <lawrence.mitchell@imperial.ac.uk> wrote:

-----BEGIN PGP SIGNED MESSAGE----- Hash: SHA1

Hi Justin,

On 14/09/15 10:57, Justin Chang wrote:

...

Hi everyone,

Attached is a code in which I am attempting to do optimization.

For the dffusion operator, I want to enforce this:

min 1/2*x^T*H*x - x^T*f s.t. 0 < x < 1

Starting at line 197 is where I have set up TAO through petsc4py. I have done this problem previously purely through PETSc (Matt's FEM w/DMPlex), and directly translated my C code into what I currently have, but the issue I am running into is enforcing non-homogeneous boundary conditions. I am still a little confused on how to enforce these in the context of firedrake/petsc4py, so I was wondering if you guys could help me "fix" what I have.

The attached code runs, but I do not think the boundary conditions for the optimization problem are not enforced correctly.

If I understand correctly, you want the resulting gradient function to satisfy the boundary conditions bc_D, yes?

If so, I think what you want to write is:

def ObjGrad_D(tao, petsc_x, petsc_g): # gradient vector - with bcs with b_D_test.dat.vec as b_D_vec: A_D_test.M.handle.mult(petsc_x,petsc_g) petsc_g.axpy(-1.0,b_D_vec) arr = petsc_g.array for bc in bc_D: # Only select nodes that are local to this process. nodes = bc.nodes[bc.nodes < bc.function_space().dof_dset.size] if isinstance(bc.function_arg, Function): # Copy over appropriate values with bc.function_arg.dat.vec_ro as bc_vec: arr[nodes] = bc_vec.array_r[nodes] else: # Scalar value, splat to all nodes arr[nodes] = float(bc.function_arg) # objective function - no bcs with b_D_nobc.dat.vec as b_D_vec: A_D_nobc.M.handle.mult(petsc_x,petsc_g) xtHx = petsc_x.dot(petsc_g) xtf = petsc_x.dot(b_D_vec) petsc_g.axpy(-1.0,b_D_vec) return 0.5*xtHx - xtf

Note I've done two things. I've swapped round the order in which you compute the gradient and the objective. Previously, you computed the objective first, but the return statement meant you never reached the gradient computation code.

The second part is to add the boundary condition application.

I loop over the boundary conditions and pull out the nodes it applies to. Then I set the value in the result array.

If you need to set homogenized bcs, then just replace the RHS assignment with just:

arr[nodes] = 0.0

Hope this helps!

Lawrence -----BEGIN PGP SIGNATURE----- Version: GnuPG v1

iQEcBAEBAgAGBQJV9qldAAoJECOc1kQ8PEYvxqwH/3Knre6BqdCqjbbyEREhYxLF TFL4ZS+tTKPqcsE8gCtw/eADKTrZk1XVaZevrWKjewLuITRNXBte2gQijlS7TggQ C+PQlhw3UTRr5yA94POmRe/bc5Ed0SUSC606WX3AI/NETp3ZnHTy/V1TBFXq5zpr hBaoJFgyJF24y+ds7TaSsuzGPmrE6yE5xEEjoqfP931wxJQ/D8N2Fxt74ohdJCN0 wiHApcE805pYieqi1SwDoLRFDBzMP83EsJ4D8B4yDnoZPW2Ow1pmL9GzRxa0iIv/ 2UbTMydpZtS+2ORMg8dG8hQGtYOkYU4YHE05j48QKUjcY7e3a533RWPp7vatnJI= =hqHs -----END PGP SIGNATURE-----

_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake

Lawrence, I want to replicate what I have in this paper in Firedrake http://arxiv.org/pdf/1506.08435.pdf In that paper I used Matt's PetscFE DMPlex directly. In his framework he removes the constrained dofs, so the parallel solvers only work with the free dofs. One of the core principles to our work is that the boundary conditions always satisfy the discrete maximum principles, that is they *are* the upper and lower bounds that the free dofs may potentially violate. So when we optimize the problem, we only optimize the free dofs. As for the way you guys do boundary conditions, I would think just applying bcs to H and f would suffice since it's basically the same as solving Hx = f. I could be wrong though, may have to do a little research first. On Tuesday, September 15, 2015, Lawrence Mitchell < lawrence.mitchell@imperial.ac.uk> wrote:

-----BEGIN PGP SIGNED MESSAGE----- Hash: SHA1

Hi Justin,

On 15/09/15 09:49, Justin Chang wrote:

...

Lawrence,

Thanks for all of your help. Here’s what I ended up with: —————— def ObjGrad_D(tao, petsc_x, petsc_g): with b_D_nobc.dat.vec as b_D_vec: A_D_nobc.M.handle.mult(petsc_x,petsc_g) xtHx = petsc_x.dot(petsc_g) xtf = petsc_x.dot(b_D_vec) petsc_g.axpy(-1.0,b_D_vec) arr = petsc_g.array for bc in bc_D: nodes = bc.nodes[bc.nodes < bc.function_space().dof_dset.size] if isinstance(bc.function_arg, Function): with bc.function_arg.dat.vec_ro as bc_vec: arr[nodes] = bc_vec.array_r[nodes] else: arr[nodes] = float(bc.function_arg) return 0.5*xtHx - xtf ——————— However it seems something is still off about the boundary conditions I am applying. I have attached three different scenarios for my two solutions psi_A and psi_B. The boundary conditions are:

psi_A: set to 1 at x == 0.0 and 0.0 <= y <= 0.5, 0 everywhere else psi_B: set to 1 at x == 0.0 and 0.5 < y <= 1.0, 0 everywhere else

- no_opt_psiA/B.png illustrates solving the standard diffusion without any optimization. As far as I know, these solutions are “correct” in terms of the boundary conditions.

- opt_wo_bounds_psiA/B.png illustrates solving the standard diffusion with TAO’s optimization but without any variable bounds. I would expect this to have the same solution as the no_opt cases but it seems to me the boundary conditions are messed up; face ID’s 3 and 4 (i.e., top and bottom) should have homogeneous bcs but it seems they are not registering.

- opt_w_bounds_psiA/B.png now includes the variable bounds set by these lines:

lb = Function(W) ub = Function(W) lb.assign(0) ub.assign(1) tao = PETSc.TAO().create(MPI.COMM_WORLD) with lb.dat.vec as lb_vec, ub.dat.vec as ub_vec: tao.setVariableBounds(lb_vec,ub_vec)

The boundary conditions are also screwed up at the top and bottom faces. In this case, psi_B (corresponding to W.sub(1)) is even more screwed up - the non-homogeneous part of psi_B isn’t even registering (hence concentrations of < 1).

Replacing the RHS of the arr[nodes] assignments with 0 didn’t change the solution either. Any idea what’s going on here?

I don't have any immediate insights, sorry. But some comments. When you're minimising

0.5 x^T H x - x^T f

What parts of that do you expect to obey the boundary conditions? I presume that f should? I'm a bit lost in all the code. Can you maybe write down the maths (with the specification of the bcs) and we can go from them instead?

Lawrence -----BEGIN PGP SIGNATURE----- Version: GnuPG v1

iQEcBAEBAgAGBQJV9+zQAAoJECOc1kQ8PEYvex8IAI+biiyYkJmhE57nyS13tiNs PxS/TLu1UFMAHlXs272S3r+U4DMXeVFHeb2kAf9hwSUf9g25+38rAfSYX2SH5KX2 +uXpCeOJAuxquS/zI9i7y8AP+ftbTMbmJOGLjnFZRlEDsHVjvqEToPzO13JnXCYp SIpGId6aO0FbFwT94ZDIV9fQKn2bYsLw9FSjWAEVCXoZiL6yEeyNstH32jSf8EU4 JrXsyekL4IGc8t5UtfpWiumVRMwRdfQntA+oqpEyeKiqDLFdTACAMzuAxDyXwbXH GFLcN3yYMU52R82KpdsymOEodYRaliF4CvULXGmEQwJkmoXJJvDUFnPF0b37VJ0= =JH/J -----END PGP SIGNATURE-----

_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk <javascript:;> https://mailman.ic.ac.uk/mailman/listinfo/firedrake

Okay I just figured it out. When assembling the A_D and b_D terms (corresponding to H and f in the convex optimization problem), only apply bcs to b_D. Originally I was applying bc_D to both A_D and b_D, but I guess the solution was to apply it to b_D only. Anyway, thanks again for your help :) Justin On Tue, Sep 15, 2015 at 4:39 AM, Justin Chang <jychang48@gmail.com> wrote:

Also,

The code I have right is similar to the problem described in Section 5.2 of this paper

http://arxiv.org/pdf/1210.5290.pdf

I am not concerned with the reactions at the moment, only the non negativity of the diffusion operator.

Thanks, Justin

On Tuesday, September 15, 2015, Justin Chang <jychang48@gmail.com> wrote:

...

Lawrence,

I want to replicate what I have in this paper in Firedrake

http://arxiv.org/pdf/1506.08435.pdf

In that paper I used Matt's PetscFE DMPlex directly. In his framework he removes the constrained dofs, so the parallel solvers only work with the free dofs. One of the core principles to our work is that the boundary conditions always satisfy the discrete maximum principles, that is they *are* the upper and lower bounds that the free dofs may potentially violate. So when we optimize the problem, we only optimize the free dofs. As for the way you guys do boundary conditions, I would think just applying bcs to H and f would suffice since it's basically the same as solving Hx = f. I could be wrong though, may have to do a little research first.

On Tuesday, September 15, 2015, Lawrence Mitchell < lawrence.mitchell@imperial.ac.uk> wrote:

...

-----BEGIN PGP SIGNED MESSAGE----- Hash: SHA1

Hi Justin,

On 15/09/15 09:49, Justin Chang wrote:

...

Lawrence,

Thanks for all of your help. Here’s what I ended up with: —————— def ObjGrad_D(tao, petsc_x, petsc_g): with b_D_nobc.dat.vec as b_D_vec: A_D_nobc.M.handle.mult(petsc_x,petsc_g) xtHx = petsc_x.dot(petsc_g) xtf = petsc_x.dot(b_D_vec) petsc_g.axpy(-1.0,b_D_vec) arr = petsc_g.array for bc in bc_D: nodes = bc.nodes[bc.nodes < bc.function_space().dof_dset.size] if isinstance(bc.function_arg, Function): with bc.function_arg.dat.vec_ro as bc_vec: arr[nodes] = bc_vec.array_r[nodes] else: arr[nodes] = float(bc.function_arg) return 0.5*xtHx - xtf ——————— However it seems something is still off about the boundary conditions I am applying. I have attached three different scenarios for my two solutions psi_A and psi_B. The boundary conditions are:

psi_A: set to 1 at x == 0.0 and 0.0 <= y <= 0.5, 0 everywhere else psi_B: set to 1 at x == 0.0 and 0.5 < y <= 1.0, 0 everywhere else

- no_opt_psiA/B.png illustrates solving the standard diffusion without any optimization. As far as I know, these solutions are “correct” in terms of the boundary conditions.

- opt_wo_bounds_psiA/B.png illustrates solving the standard diffusion with TAO’s optimization but without any variable bounds. I would expect this to have the same solution as the no_opt cases but it seems to me the boundary conditions are messed up; face ID’s 3 and 4 (i.e., top and bottom) should have homogeneous bcs but it seems they are not registering.

- opt_w_bounds_psiA/B.png now includes the variable bounds set by these lines:

lb = Function(W) ub = Function(W) lb.assign(0) ub.assign(1) tao = PETSc.TAO().create(MPI.COMM_WORLD) with lb.dat.vec as lb_vec, ub.dat.vec as ub_vec: tao.setVariableBounds(lb_vec,ub_vec)

The boundary conditions are also screwed up at the top and bottom faces. In this case, psi_B (corresponding to W.sub(1)) is even more screwed up - the non-homogeneous part of psi_B isn’t even registering (hence concentrations of < 1).

Replacing the RHS of the arr[nodes] assignments with 0 didn’t change the solution either. Any idea what’s going on here?

I don't have any immediate insights, sorry. But some comments. When you're minimising

0.5 x^T H x - x^T f

What parts of that do you expect to obey the boundary conditions? I presume that f should? I'm a bit lost in all the code. Can you maybe write down the maths (with the specification of the bcs) and we can go from them instead?

Lawrence -----BEGIN PGP SIGNATURE----- Version: GnuPG v1

iQEcBAEBAgAGBQJV9+zQAAoJECOc1kQ8PEYvex8IAI+biiyYkJmhE57nyS13tiNs PxS/TLu1UFMAHlXs272S3r+U4DMXeVFHeb2kAf9hwSUf9g25+38rAfSYX2SH5KX2 +uXpCeOJAuxquS/zI9i7y8AP+ftbTMbmJOGLjnFZRlEDsHVjvqEToPzO13JnXCYp SIpGId6aO0FbFwT94ZDIV9fQKn2bYsLw9FSjWAEVCXoZiL6yEeyNstH32jSf8EU4 JrXsyekL4IGc8t5UtfpWiumVRMwRdfQntA+oqpEyeKiqDLFdTACAMzuAxDyXwbXH GFLcN3yYMU52R82KpdsymOEodYRaliF4CvULXGmEQwJkmoXJJvDUFnPF0b37VJ0= =JH/J -----END PGP SIGNATURE-----

_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake

-----BEGIN PGP SIGNED MESSAGE----- Hash: SHA1 On 21/09/15 10:38, Justin Chang wrote:

Hello again,

I am still having issues with this, what I had done originally was not correct.

In FEniCS, all I had to do was

A_D, b_D = assemble_system(a_D, L_D, bcs_D)

And feed the A_D and b_D terms into scipy's optimization. Essentially I am solving:

min 0.5*x^T*A_D*x - x^T*b_D

and the gradient of the above would be:

A_D*x - b_D

I have attached said FEniCS code (named "Decay_Problem_3.py") which solves a different problem, but the concept is the same.

However, this methodology does not work with my current setup. If I replicate the above in Firedrake, the entire solution seems to be fixated to the dirichlet BCs and whatever initial conditions I have provided, which I am guessing is what you were warning me about previously. Do you guys handle Dirichlet boundary conditions differently than how FEniCS handles theirs? I have also attached my current firedrake code (named "2D_plume_ADR_ex1.py", the parts in question begin at line 161)

Ah, OK. So when you assemble an operator with boundary conditions, the bc rows /and/ columns are zeroed and a 1 is placed on the diagonal. The appropriate RHS for this linear is system is then not just the linear form L_D assembled with the bcs applied. You need to "lift" the boundary conditions from the operator to the RHS. assemble_system does this under the hood for you. I think things will work if you do the following: # I assume the solution space is W A_D = assemble(a_D, bcs=bcs_D) tmp = Function(W) for bc in bcs_D: bc.apply(tmp) rhs_bcs = assemble(action(a_D, tmp)) b_D = assemble(L_D) b_D -= rhs_bcs for bc in bcs_D: bc.apply(b_D) More details on the background here: www.firedrakeproject.org/boundary_conditions.html Cheers, Lawrence -----BEGIN PGP SIGNATURE----- Version: GnuPG v1 iQEcBAEBAgAGBQJV/9K0AAoJECOc1kQ8PEYv0e0H/0OPFWJlgURq0MDdW6322f3I 3lfPOZH8Odlj4sd4XX1+bx9/q2dN9Ob1reWFFWRalVukIS6JsJMG4All07X1ZZwi qyV3j4fZNoucmzf5taaM46EoYNE+IS7COKrfwhMU+jS+erPYL0JdGsZHneUazjet /wt3MrMKt5Njlzw39yVKAMpD7WOXda/568lcEWpdce3nR+tEjhdK92v9DlSt1CdQ /GoiSgmAAc5PHguOqeGcbcXasp4cqvoiyMv5Mehk5gIdwO747c0jruoxQK8Mx3V8 uBkfalDAUOhgMKgllPESD3gFAJbjopx2MbB0W8bEPYz2n4tgb/ZOujqb6ScZfrY= =bVZX -----END PGP SIGNATURE-----