Subject: Re: [Ffc] [Branch ~ffc-core/ffc/main] Rev 1369: Fixed Garth's optimised quadrature problem.

From: Anders Logg <[email-masked]>
Date: Mon, 1 Feb 2010 21:37:51 +0100

Impressive! :-)

Is it release time or are we waiting for Kristian to do some magic
with evaluate_derivatives for Hdiv/curl?

--
Anders

On Mon, Feb 01, 2010 at 08:30:28PM -0000, [email-masked] wrote:
> ------------------------------------------------------------
> revno: 1369
> committer: Marie E. Rognes <[email-masked]>
> branch nick: ffc
> timestamp: Mon 2010-02-01 21:27:04 +0100
> message:
>   Fixed Garth's optimised quadrature problem.
> modified:
>   ffc/quadrature/fraction.py
>
>

> === modified file 'ffc/quadrature/fraction.py'
> --- ffc/quadrature/fraction.py	2010-01-30 19:50:03 +0000
> +++ ffc/quadrature/fraction.py	2010-02-01 20:27:04 +0000
> @@ -5,7 +5,7 @@
>  __copyright__ = "Copyright (C) 2009-2010 Kristian B. Oelgaard"
>  __license__  = "GNU GPL version 3 or any later version"
>
> -# Last changed: 2010-01-21
> +# Last changed: 2010-02-01
>
>  # FFC modules.
>  from ffc.log import error
> @@ -241,9 +241,14 @@
>          self = found*remain."""
>
>          # Reduce the numerator by the var type.
> -        # Handle case where numerator is a sum.
>          if self.num._prec == 3:
> -            num_found, num_remain = self.num.reduce_vartype(var_type)[0]
> +            foo = self.num.reduce_vartype(var_type)
> +            if len(foo) == 1:
> +                num_found, num_remain = self.num.reduce_vartype(var_type)[0]
> +            else:
> +                # meg: I have only a marginal idea of what I'm doing here!
> +                return create_sum([create_fraction(create_product([num_found, num_remain]), self.denom)
> +                                   for (num_found, num_remain) in foo]).expand().reduce_vartype(var_type)
>          else:
>              num_found, num_remain = self.num.reduce_vartype(var_type)
>
> @@ -259,7 +264,7 @@
>              denom_found, denom_remain = self.denom.reduce_vartype(var_type)
>
>          # If we have a Sum in the denominator, all terms must be reduced by
> -        # the same terms to make sense.
> +        # the same terms to make sense
>          else:
>              remain = []
>              for m in self.denom.vrs:
> @@ -294,6 +299,7 @@
>          # TODO: Add more checks to avoid expansion.
>          found = None
>          # There is always a remainder.
> +
>          remain = create_fraction(num_remain, denom_remain).expand()
>
>          if num_found:
> @@ -306,6 +312,7 @@
>                  found = create_fraction(create_float(1), denom_found)
>              else:
>                  found = ()
> +
>          return (found, remain)
>
>  # FFC quadrature modules.
>

---

From: Marie Rognes <[email-masked]>
Date: Mon, 01 Feb 2010 21:39:07 +0100

> with evaluate_derivatives for Hdiv/curl?

Who the heck needs evaluate_derivatives for Hdiv/curl ;) ?

Throw error?

--
Marie

> ------------------------------------------------------------------------
>
> _______________________________________________
> Mailing list: https://launchpad.net/~ffc
> Post to     : [email-masked]
> Unsubscribe : https://launchpad.net/~ffc
> More help   : https://help.launchpad.net/ListHelp
>

---

From: Anders Logg <[email-masked]>
Date: Mon, 1 Feb 2010 21:41:23 +0100

> Who the heck needs evaluate_derivatives for Hdiv/curl ;) ?

Yes, if it looks like it will take time to fix. It hasn't worked
before as far as I understand so we don't necessarily need to make
it work now. Kristian?

--
Anders

---

From: Kristian Oelgaard <[email-masked]>
Date: Mon, 1 Feb 2010 21:49:59 +0100 (CET)

> Impressive! :-)

It sure is, many things could be handled in a much simpler way if it wasn't for memory and runtime considerations of the optimisations. If I get around to implementing the other optimisation strategy, we could make life simpler in these classes and only let the symbolics stuff work on moderately complex forms.

> with evaluate_derivatives for Hdiv/curl?

Actually, I'm waiting for Marie to fix this one too :)
I'm having problems with the transformation back to the physical element, at least that's where I think the problem is. On the other hand, when I'm computing values on the reference element (for the tests) should it really matter for these elements? I don't see any transformation in FIAT for Hdiv/curl elements.

Kristian

> -----BEGIN PGP SIGNATURE-----
> Version: GnuPG v1.4.9 (GNU/Linux)
>
> iEYEARECAAYFAktnO54ACgkQTuwUCDsYZdEMIgCfQj7hbvwXKDVzP72ykidr29E9
> 3s8An1xmHf6nfXJUbnwkwCeZ7dDq5cO0
> =9wkW
> -----END PGP SIGNATURE-----
>

---

From: Marie Rognes <[email-masked]>
Date: Mon, 01 Feb 2010 22:02:05 +0100

Kristian Oelgaard wrote:
> It sure is, many things could be handled in a much simpler way if it 
> wasn't for memory and runtime considerations of the optimisations. If 
> I get around to implementing the other optimisation strategy, we could 
> make life simpler in these classes and only let the symbolics stuff 
> work on moderately complex forms.
>

Sorry if I messed something up. I just thought it was better to have 
something working (but suboptimal) than something broken (and optimal ;) )

> Actually, I'm waiting for Marie to fix this one too :)

I'm taking a look, but it is not likely to happen within the next couple 
of hours.

I'm lost in the vectors: What is 'derivatives' and 'values' supposed to 
be holding? (Math, please ;))

--
Marie

> I'm having problems with the transformation back to the physical 
> element, at least that's where I think the problem is. On the other 
> hand, when I'm computing values on the reference element (for the 
> tests) should it really matter for these elements? I don't see any 
> transformation in FIAT for Hdiv/curl elements.
>
> Kristian
>
>> -- 
>> Anders
>>
>>
>>> +                return 
>>> create_sum([create_fraction(create_product([num_found, num_remain]), 
>>> self.denom)

---

From: Kristian Oelgaard <[email-masked]>
Date: Mon, 1 Feb 2010 22:23:53 +0100 (CET)

On 1 February 2010 22:02, Marie Rognes <[email-masked]> wrote:
>> wasn't for memory and runtime considerations of the optimisations. If I get
>> around to implementing the other optimisation strategy, we could make life
>> simpler in these classes and only let the symbolics stuff work on moderately
>> complex forms.
>>
>

Agree, the code looks OK.

> I'm taking a look, but it is not likely to happen within the next couple of
> hours.
>
> I'm lost in the vectors: What is 'derivatives' and 'values' supposed to be
> holding? (Math, please ;)

'values' is the argument to evaluate_basis_derivatives that will hold the return values (derivatives on the physical cell) from the function. It is the users responsibility to make sure it has the correct size. 

'derivatives' is the derivatives on the reference element, which is computed like it is done inside FIAT (multiplication of coefficients, dmats and basis values of the expansion).

Kristian
>
> --
> Marie
>
>
>> I'm having problems with the transformation back to the physical element,
>> at least that's where I think the problem is. On the other hand, when I'm
>> computing values on the reference element (for the tests) should it really
>> matter for these elements? I don't see any transformation in FIAT for
>> Hdiv/curl elements.
>>
>> Kristian
>>
>>> --
>>> Anders
>>>
>>>

---

From: Marie Rognes <[email-masked]>
Date: Mon, 01 Feb 2010 22:32:01 +0100

Kristian Oelgaard wrote:
> 'values' is the argument to evaluate_basis_derivatives that will hold 
> the return values (derivatives on the physical cell) from the 
> function. It is the users responsibility to make sure it has the 
> correct size.
> 'derivatives' is the derivatives on the reference element, which is 
> computed like it is done inside FIAT (multiplication of coefficients, 
> dmats and basis values of the expansion).
>

Ok. And the values are ordered so that first all derivatives for value 
dimension 0, then all derivatives for value dimension 1, etc.?

Currently only the first two (out of four?) 'derivatives' are mapped for 
BDM _1on triangle?

        derivatives[0] = (1.0/detJ)*((J_00*tmp_ref0 + J_01*tmp_ref1));
        derivatives[1] = (1.0/detJ)*((J_10*tmp_ref0 + J_11*tmp_ref1));

but you use higher indices later:

        values[row] += transform[row][col]*derivatives[col];
        values[num_derivatives + row] += 
transform[row][col]*(derivatives[num_derivatives + col]);

Is this intentional?

--
Marie

>>> element,
>>>>>         # If we have a Sum in the denominator, all terms must be 
>>>>> reduced
>>>>> by

---

From: Kristian Oelgaard <[email-masked]>
Date: Mon, 1 Feb 2010 22:37:39 +0100 (CET)

On 1 February 2010 22:32, Marie Rognes <[email-masked]> wrote:
>>> of
>>> be
>> 'values' is the argument to evaluate_basis_derivatives that will hold the
>> return values (derivatives on the physical cell) from the function. It is
>> the users responsibility to make sure it has the correct size.
>> computed like it is done inside FIAT (multiplication of coefficients, dmats
>> and basis values of the expansion).
>>
>

Yes. There is a comment at the top of evaluatebasisderivatives.py, but it might be a good idea to put in the generated code too?

> Currently only the first two (out of four?) 'derivatives' are mapped for BDM
> _1on triangle?
>

I pushed a fix for this a while back, I was hoping that was the error, but I still get wrong results.
It should be like below of course.

Kristian

---

From: Marie Rognes <[email-masked]>
Date: Mon, 01 Feb 2010 22:57:15 +0100

Kristian Oelgaard wrote:
>>>>> wasn't for memory and runtime considerations of the optimisations. 
>>>>> If I
>>>>> get
>>>> something working (but suboptimal) than something broken (and 
>>>> optimal ;)
>>>> )
>>>
> Yes. There is a comment at the top of evaluatebasisderivatives.py, but 
> it might be a good idea to put in the generated code too?
>

Sorry, my bad -- something strange happened.

And I guess you have checked the indices a couple of hundred times?

How about the famous UFC/FIAT reference element difference?

--
Marie

>>>>>>> +                # meg: I have only a marginal idea of what I'm 
>>>>>>> doing
>>>>>>> here!
>>>>>>> +                return

---

From: Kristian Oelgaard <[email-masked]>
Date: Mon, 1 Feb 2010 23:04:44 +0100 (CET)

On 1 February 2010 22:57, Marie Rognes <[email-masked]> wrote:
>>>>>> I
>>>>> something working (but suboptimal) than something broken (and optimal
>>>>> ;)
>>>>> )
>>>>
>>>> the
>> Yes. There is a comment at the top of evaluatebasisderivatives.py, but it
>> might be a good idea to put in the generated code too?
>>
>>> BDM

I think I found it, I was mapping values from FIAT incorrectly. Stupid, stupid me.
I'll just clean things up then Anders can make a release.

Kristian

---

From: Marie Rognes <[email-masked]>
Date: Mon, 01 Feb 2010 23:07:20 +0100

Kristian Oelgaard wrote:
>>>>>>> around to implementing the other optimisation strategy, we could 
>>>>>>> make
>>>>>>> life
>>>>>> optimal
>>>>> return values (derivatives on the physical cell) from the 
>>>>> function. It
>>>>> is
>>> I pushed a fix for this a while back, I was hoping that was the 
>>> error, but
>>> I still get wrong results.
>>
>>
> I think I found it, I was mapping values from FIAT incorrectly. 

Hurray! *crossing fingers*

--
Marie

> Stupid, stupid me.
>>>>>>> at least that's where I think the problem is. On the other hand, 
>>>>>>> when
>>>>>>> I'm
>>>>>>>>> +                                   for (num_found, 
>>>>>>>>> num_remain) in
>>>>>>>>> foo]).expand().reduce_vartype(var_type)
>>>>>>>>>       else:
>>>>>>>>>           num_found, num_remain =

---

From: Anders Logg <[email-masked]>
Date: Mon, 1 Feb 2010 23:19:26 +0100

> I'll just clean things up then Anders can make a release.

Yes! Yes! Yes! :-)

I'll start with FIAT, UFC, UFL and FFC. Then recompile all forms in
DOLFIN and wait until tomorrow with the release of DOLFIN.

We need to hurry up before Garth finds a new bug. He seems determined
to delay the release. ;-)

--
Anders

---

From: Mehdi Nikbakht <[email-masked]>
Date: Mon, 01 Feb 2010 23:19:32 +0100

> Hurray! *crossing fingers*

Congratulation for such a nice job!

Mehdi

---

From: Kristian Oelgaard <[email-masked]>
Date: Mon, 1 Feb 2010 23:22:29 +0100 (CET)

> to delay the release. ;-)

Good idea, but remember he has an extra hour compared to us so there might be a long list tomorrow.

Kristian

> iEYEARECAAYFAktnU24ACgkQTuwUCDsYZdERIgCcDUj9Yb+YyhP9ZCo1eyr8UbJ7
> Rt0An1/1tQmWYIyh9Yht4SgAQMRAHZ1r
> =tgJH

---

From: Anders Logg <[email-masked]>
Date: Mon, 1 Feb 2010 23:25:18 +0100

> Good idea, but remember he has an extra hour compared to us so there might be a long list tomorrow.

That's what I'm afraid of. But it also means he has an extra hour to
fix those bugs. :-)

--
Anders

---

From: "Garth N. Wells" <[email-masked]>
Date: Tue, 02 Feb 2010 00:07:08 +0000

> to delay the release. ;-)

You're lucky that I've been busy at the college this evening. Gown +
laptop would be frowned upon.

Garth