Subject: [FFC-dev] bug in simplify?

From: Kristian Oelgaard <[email-masked]>
Date: Tue, 18 Mar 2008 17:04:30 +0100

Hi,

I think there is a bug in simplify.py. Consider the following simple forms and
their output:

element = VectorElement("Discontinuous Lagrange", "triangle", 0, 4)
v = TestFunction(element)
u = TrialFunction(element)

def Mat(v):
  return [[v[0], v[1]],
          [v[2], v[3]]]

def Sym(v):
  return 0.5*(Mat(v) + transp(Mat(v)))

a = dot(Mat(v), Sym(u))*dx

with simplify:

    // Compute geometry tensors
    const double G0_ = det*(1.0 + 1.0);
    const double G1_ = det;
    const double G2_ = det;
    
    // Compute element tensor
    A[0] = 0.25*G0_;
    A[1] = 0;
    A[2] = 0;
    A[3] = 0;
    A[4] = 0;
    A[5] = 0.25*G0_;
    A[6] = 0.25*G1_;
    A[7] = 0;
    A[8] = 0;
    A[9] = 0.25*G2_;
    A[10] = 0.25*G0_;
    A[11] = 0;
    A[12] = 0;
    A[13] = 0;
    A[14] = 0;
    A[15] = 0.25*G0_;

or [2 0 0 0]
   [0 2 1 0]
   [0 1 2 0]
   [0 0 0 2]

without simplify:

    // Compute geometry tensors
    const double G0_ = det*(1.0 + 1.0);
    const double G1_ = det;
    const double G2_ = det;
    const double G3_ = det;
    const double G4_ = det;
    const double G5_ = det*(1.0 + 1.0);
    
    // Compute element tensor
    A[0] = 0.25*G0_;
    A[1] = 0;
    A[2] = 0;
    A[3] = 0;
    A[4] = 0;
    A[5] = 0.25*G1_;
    A[6] = 0.25*G2_;
    A[7] = 0;
    A[8] = 0;
    A[9] = 0.25*G4_;
    A[10] = 0.25*G3_;
    A[11] = 0;
    A[12] = 0;
    A[13] = 0;
    A[14] = 0;
    A[15] = 0.25*G5_;

or [2 0 0 0]
   [0 1 1 0]
   [0 1 1 0]
   [0 0 0 2]

which I what one would expect. So I think this indicates a bug somewhere...

Kristian

---

From: Marie Rognes <[email-masked]>
Date: Tue, 18 Mar 2008 12:06:05 -0400

> I think there is a bug in simplify.py. Consider the following simple forms and

What? A bug in simplify? Impossible ;)

I'll look at it a bit later this afternoon.

-- 
Marie E. Rognes
Ph.D Fellow, 
Centre of Mathematics for Applications, 
University of Oslo
http://folk.uio.no/meg

---

From: Kristian Oelgaard <[email-masked]>
Date: Tue, 18 Mar 2008 17:17:41 +0100

> What? A bug in simplify? Impossible ;)

I know, I almost couldn't believe it myself :)
 
> I'll look at it a bit later this afternoon.

OK

Kristian

---

From: Marie Rognes <[email-masked]>
Date: Tue, 18 Mar 2008 13:09:55 -0400

Kristian Oelgaard wrote:
>>> I think there is a bug in simplify.py. Consider the following simple forms
>>>       
>> and
>>     
>>> their output:
>>>
>>>   
>>>       

However, this time it is actually not simplify that is the main culprit:

    After index reassignment:
     
          0.5 | vi0[0]*vi1[0]*dX(0) + 0.5 | vi0[0]*vi1[0]*dX(0) + 0.5 |
    vi0[1]*vi1[1]*dX(0) + 0.5 | vi0[1]*vi1[2]*dX(0) + 0.5 |
    vi0[2]*vi1[2]*dX(0) + 0.5 | vi0[2]*vi1[1]*dX(0) + 0.5 |
    vi0[3]*vi1[3]*dX(0) + 0.5 | vi0[3]*vi1[3]*dX(0)
     
      After simplification:
     
          0.5 | vi0[b0]*vi1[b0]*dX(0) + 0.5 | vi0[b0]*vi1[b0]*dX(0) +
    0.5 | vi0[1]*vi1[2]*dX(0) + 0.5 | vi0[2]*vi1[1]*dX(0)

where the index b0 has range [0, 1, 2, 3] in the first term and [0, 3] 
in the second term. This makes sense. However, in the factorization step 
for the tensor representation, these two terms are considered equal 
since their signatures are the same and thus a faulty tensor 
representation is generated.

This is an inherent complication with regard to letting indices be 
interpreted as equivalent for the sake of index contractions vs 
different for factorizations.

I am not sure of what the best way of approaching this would be. I don't 
feel like messing too much with the factorization step in the 
tensorrepresentation...

However, I was planning on refactoring the simplify code at some point. 
I think a key part of this would be to generalize the signatures so that 
we could use these instead of elaborate checks that are there now... 
Will not be before spring break though.

--
Marie

---

From: Kristian Oelgaard <[email-masked]>
Date: Tue, 18 Mar 2008 18:21:08 +0100

Quoting Marie Rognes <[email-masked]>:

> in the second term. This makes sense.

Maybe, on a side note: the fact that 'b0' can have a different range in each
monomial makes it pretty difficult to debug/understand what is going on by just
looking at the pretty print version of the form. The index should be different
in my opinion.

 However, in the factorization step 
> feel like messing too much with the factorization step in the 

No, me neither.

> Will not be before spring break though.

Sounds like a good idea, maybe we should switch off simplify until then?

Kristian

---

From: Marie Rognes <[email-masked]>
Date: Tue, 18 Mar 2008 14:37:49 -0400

> looking at the pretty print version of the form. The index should be different

Agree. Maybe we could have two versions of printing indices, one with 
the ranges (if set) ala

    b0 (0, 1, 2, 3)

for debug purposes and a denser one for non-debug purposes?

Actually, that would be a starting point for an accurate factorization 
check too.

---

From: Anders Logg <[email-masked]>
Date: Wed, 19 Mar 2008 09:50:36 +0100

> Actually, that would be a starting point for an accurate factorization 

See if it works now. It was a simple fix (assuming it works) to
include the index range in the computation of signatures (see
signature.py).

-- 
Anders

---

From: Kristian Oelgaard <[email-masked]>
Date: Wed, 19 Mar 2008 12:08:16 +0100

> include the index range in the computation of signatures (see

Yes, it works, but the output could be optimised:

    const double G0_ = det;
    const double G1_ = det;
    const double G2_ = det;
    const double G3_ = det;
    
    // Compute element tensor
    A[0] = 0.25*G0_ + 0.25*G1_;
    A[1] = 0;
    A[2] = 0;
    A[3] = 0;
    A[4] = 0;
    A[5] = 0.25*G0_;
    A[6] = 0.25*G2_;
    A[7] = 0;
    A[8] = 0;
    A[9] = 0.25*G3_;
    A[10] = 0.25*G0_;
    A[11] = 0;
    A[12] = 0;
    A[13] = 0;
    A[14] = 0;
    A[15] = 0.25*G0_ + 0.25*G1_;

So something in factorization should change.

And it would be nice to include some information about the range of the indices
in the form signature:

  After simplification:

      0.5 | vi0[b0]*vi1[b0]*dX(0) + 0.5 | vi0[b0]*vi1[b0]*dX(0) + 0.5 |
vi0[1]*vi1[2]*dX(0) + 0.5 | vi0[2]*vi1[1]*dX(0)

b0 in first product 0-3, in second 0-1.

Kristian

> 
> -- 
> Anders
> _______________________________________________
> FFC-dev mailing list
> [email-masked]
> http://www.fenics.org/mailman/listinfo/ffc-dev
>

---

From: Anders Logg <[email-masked]>
Date: Tue, 25 Mar 2008 00:15:13 +0100

> So something in factorization should change.

Yes, but I don't want to mess with it unless I really have to.
A lot of work needs to be done when we move to UFL anyway, so
I'd rather wait.

> b0 in first product 0-3, in second 0-1.

ok, see if you like it now.

-- 
Anders

---

From: Kristian Oelgaard <[email-masked]>
Date: Tue, 25 Mar 2008 10:49:34 +0100

Quoting Anders Logg <[email-masked]>:

> > > > > Maybe, on a side note: the fact that 'b0' can have a different range
> in
> > > each
> > > > > monomial makes it pretty difficult to debug/understand what is going
> on
> > > by just

OK.
 
> ok, see if you like it now.

The output is definitely much more informative without being cluttered.

Kristian