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