[firedrake] Enforcing slip boundary condition strongly

[David Ham]

That's quite yuck, however I can see the attraction as we're kind of sort
of not mucking with the PyOP2 interface.

However the generated code is specific to BCs applied to a particular
component, which doesn't satisfy Justin's (very standard) use case.

I think we need to encode *which* component we're changing into the map
value. EG.

val = -val -(2**(30-i))

where i is the component. This can be decoded to zero the correct entry in
a general purpose way. The downside is that we limit ourselves to about
half a billion nodes per core. I think that's OK.


On Fri, 10 Jul 2015 at 11:11 Mitchell, Lawrence <
lawrence.mitchell at imperial.ac.uk> wrote:

>
> > On 10 Jul 2015, at 10:25, David Ham <David.Ham at imperial.ac.uk> wrote:
> >
> > Hi Colin,
> >
> > That does not appear to be what Justin is asking for. He is asking how
> to set a Dirichlet condition on one component of a VectorFunctionSpace. I'm
> not sure we can do that right now, but I think this should be fixable. I
> want to check what the Dolfin syntax for this is (so we maintain
> compatibility) and then try to fix it. (Of course someone might want to
> pleasantly surprise me and tell me we can do this already!)
>
> We can't do it at the moment, the rest of this email explains why and
> sketches a possible solution.
>
> OK, so there are two aspects to this, which make things somewhat tricky.
>
> Our data model for indexing vector fields is that we store the Map to the
> node and stack the dofs on top.
>
> To apply such slip BCs, where we only set one of the two dofs, we need to
> do the following:
>
> - Apply the BC value to the appropriate part of the vector.  I think this
> is reasonably straightforward and appear to have something that works.
>
> - Fix up the operator appropriately.  This latter option is trickier.  The
> way we currently do this is to drop boundary conditions on the floor when
> assembling the operator.  This is done by flagging the appropriate
> rows/columns to be dropped by putting a negative value in the Map at the
> appropriate place and letting PETSc do the rest.
>
> For vector-valued fields we do block insertion, i.e. we insert a 2x2 block
> at each location in the matrix indexed by the base map.  So the current
> setup can only drop complete blocks, rather than part of one (as required
> in this case).
>
> The question is, then, how to communicate the relevant bit of information
> into the code that inserts.
>
> How about this:
>
> When we splat the values in the Map we convert:
>
> val => -(val + 1)
>
> If we're apply the BC on part of the map we decorate it with the index
> offset.
>
> We then build a local (unrolled map) and fix up the entries.  The negative
> offset encoding allows us to recover the correct value.
>
> So assembly turns from:
>
> A = build_element_tensor(...);
> MatSetValuesBlockedLocal(block_rows, block_cols, A);
>
> Into:
>
> A = build_element_tensor(...);
> {
>     PetscInt rows[arity*dim];
>     PetscInt cols[arity*dim];
>     // Assume we're setting the bc on the second component.
>     for ( i = 0; i < arity; i++ ) {
>          if (row_map[i] < 0) {
>              rows[i*dim + 0] = dim*(-(row_map[i] + 1));
>              rows[i*dim + 1] = -1;
>          } else {
>               rows[i*dim + 0] = dim*row_map[i];
>               rows[i*dim + 1] = dim*row_map[i] + 1;
>          }
>     }
>     // Same for cols
>     MatSetValuesLocal(rows, cols, A);
> }
>
> All the other options I came up with seem more complicated.
>
> Lawrence
>
-------------- next part --------------
HTML attachment scrubbed and removed