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 


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