Treat transformation coefficients (J, invJ, detJ) specially

[kynan]

It is necessary to distinguish between coefficients representing a (physical) field and coefficients representing output from transformation-to-physical (J/invJ/detJ):

• Field coefficients are defined on nodes of the mesh and hence need to be evaluated at the quadrature points (in thread local on-chip memory) before being used in the quadrature loop to evaluate the integrals.
• J/invJ/detJ are defined on the elements instead and this evaluation is not necessary. Values are retrieved from global device memory:
  • if space is linear, there is one value per element to be used at all the quadrature points
  • otherwise there is one value per quadrature point which is to be used

    Special treatment of J/invJ/detJ compared to field coefficients

In the form, J/invJ/detJ differ from field coefficients as follows:

• no evaluation at quadrature points is to be generated
• no thread-local on-chip memory is allocated
• the subscript needs to index into global memory, taking into account element id etc.

In the assembler, J/invJ/detJ differ from field coefficients as follows:

• there is a call to the transformation kernel with J/invJ/detJ as output parameters
• J/invJ/detJ do not interact with state: neither are they extracted from state nor written back to state
• global device memory needs to be allocated separately for J/invJ/detJ, since this is not done by state

  Distinction in UFL and code generator

• J/invJ/detJ could appear in the UFL as Coefficient with a special count (-3/-4/-5). However, that gets lost when coefficients are renumbered in the canonicalising step. Hence the Form and Assembler backends would need to go back to the original coefficient attached to the form_data go get the original count.
• There could subclasses Jacobian, JacobianInverse and JacobianDeterminant for J/invJ/detJ, which may make things easier. This would also allow attaching further meta data later on (which is not possible for Coefficient since it's protected by __slots__). This could be used e.g. to indicate whether the evaluation is to be done in the tensor tabulation kernel or in a separate transformation kernel.

[kynan]

After a discussion with Simon I think most of these arguments don't apply and there may be no need for any special treatment:

• Having J/invJ/detJ in state may actually be a very good idea, allowing to cache the Jacobian if it doesn't change over time. State would then also take care of device memory (de)allocations.
• For the evaluation at quadrature points, the same rules as for other coefficients apply (i.e. it is not necessary to do any evaluation if the coefficient lives on P₀^DG

[dham]

I don't think I agree with this.

J, invJ and detJ are just coefficients like any other. It is true that in the linear element case they are constant for each element but that is true for any P0 coefficient. The optimisations described below are a good idea but they can be applied simply by noticing which fields are P0.

On 17 August 2011 14:49, kynan reply@reply.github.com wrote:

It is necessary to distinguish between coefficients representing a (physical) field and coefficients representing output from transformation-to-physical (J/invJ/detJ):

• Field coefficients are defined on nodes of the mesh and hence need to be evaluated at the quadrature points (in thread local on-chip memory) before being used in the quadrature loop to evaluate the integrals.
• J/invJ/detJ are defined on the elements instead and this evaluation is not necessary. Values are retrieved from global device memory:  * if space is linear, there is one value per element to be used at all the quadrature points  * otherwise there is one value per quadrature point which is to be used

Special treatment of J/invJ/detJ compared to field coefficients

In the form, J/invJ/detJ differ from field coefficients as follows:

• no evaluation at quadrature points is to be generated
• no thread-local on-chip memory is allocated
• the subscript needs to index into global memory, taking into account element id etc.

In the assembler, J/invJ/detJ differ from field coefficients as follows:

• there is a call to the transformation kernel with J/invJ/detJ as output parameters
• J/invJ/detJ do not interact with state: neither are they extracted from state nor written back to state
• global device memory needs to be allocated separately for J/invJ/detJ, since this is not done by state

Distinction in UFL and code generator

• J/invJ/detJ could appear in the UFL as Coefficient with a special count (-3/-4/-5). However, that gets lost when coefficients are renumbered in the canonicalising step. Hence the Form and Assembler backends would need to go back to the original coefficient attached to the form_data go get the original count.
• There could subclasses Jacobian, JacobianInverse and JacobianDeterminant for J/invJ/detJ, which may make things easier. This would also allow attaching further meta data later on (which is not possible for Coefficient since it's protected by __slots__). This could be used e.g. to indicate whether the evaluation is to be done in the tensor tabulation kernel or in a separate transformation kernel.

Reply to this email directly or view it on GitHub: https://github.com/gmarkall/manycore_form_compiler/issues/38

[kynan]

Which part do you not agree with? As clarified in the last post I think no special treatment is necessary and the same optimisations can be applied as to any other coefficients.

Do you think subclassing is necessary / a good idea?

[dham]

I think there is a misunderstaning here. I didn't get the email for your post of August 17 so what I was not agreeing with is the original idea of special-casing Jacobians. I agree that they should not be treated as special.

[kynan]

This issue and the discussion assume that J/detJ/invJ are computed by a separate kernel and passed to the form evaluation kernel as parameters. Since we are computing those inside the kernel, this is not relevant.