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
