-----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-----