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