ianhinder / Kranc

A Mathematica package for generating code for solving time dependent partial differential equations
http://kranccode.org
GNU General Public License v2.0
28 stars 10 forks source link

Inconsistent indices don't lead to error #108

Open eschnett opened 10 years ago

eschnett commented 10 years ago

I had the following expression in a calculation:

    dot[B[ua]] -> IfB[+ dotXt[ua]
                      + IfThen[fixAdvectionTerms!=0,
                               Upwind[beta[ub], Xt[ua], lb],
                               0]
                      - (betaDriverValue
                         (+ B[ua]
                          + IfThen[fixAdvectionTerms!=0 && advectShift!=0,
                                   Upwind[beta[ub], beta[ua], la] /
                                   (shiftGammaCoeffValue alpha^shiftAlphaPower),
                                   0]))
                      + IfThen[fixAdvectionTerms==0 && advectShift!=0,
                               Upwind[beta[ub], B[ua], lb],
                               0]
                      + Dissipation[B[ua]],
                      0],

Note that the index "la" appears in there (in the second Upwind term), which it should not, as the respective subexpression then has index "ub" free, inconsistent with the left-hand side. "la" needs to be replaced by "lb".

This did not lead to an error.

Upwind is defined as follows:

Upwind[dir_, u_,il_] := dir PDua[u,il] + Abs[dir] PDus[u,il];

Is there a way to declare functions taking tensor indices?

ianhinder commented 10 years ago

I believe the problem is that the checking code doesn't know what to make of IfB, or indeed IfThen. As far as TensorTools is concerned, these are just functions. How should TensorTools handle unknown functions? Should it insist that each argument has the same free indices? Sometimes the first argument of IfThen will have the same indices as the second and third argument, but sometimes it will be a scalar expression. It looks like TensorTools needs extra information about specific functions in order to know what is valid.

eschnett commented 10 years ago

There are some expressions that "eat tensors", such as MatrixOfComponents. In other cases, the functions should probably be treated as element-wise non-linear transformations of tensor components. In the case of IfThen[a,b,c], the condition a should not have any tensor indices, and the branches b and c must conform to each other. In this sense, "conform" means that scalars may be promoted to tensors. Similar properties hold for Max[a,b] and Min[a,b].

barrywardell commented 10 years ago

I encountered a similar IfThen problem recently while tidying up the xTensor support in Kranc. Since xTensor does more strict error checking than TensorTools, and xTensor didn't know how to handle IfThen, it complained with some standard Kranc scripts such as McLachlan (I think it would also have picked up the problem mentioned above).

The solution I came up with was the following:

  1. "Factor out" the conditionals within the right hand side. What I mean by this is that I turned all possible branches (since there can be multiple, possibly nested IfThen's) into a single switch, with one right hand side for each possible branch of the switch.
  2. Check each individual right hand side for tensor-correctness (including the left hand side)
  3. Expand out components for each right-hand side separately.
  4. Turn the result into the form, where all possible branches are represented:
lhs -> IfThen[condition1 && condition2, rhs1true2true, IfThen[condition1 && !condition2, rhs1true2false, ...]]

Step 4 might also be better implemented in some other way, e.g. pulling the conditional out further to get a set like IfThen[condition1 && condition2, lhs -> rhs1true2true, IfThen[condition1 && !condition2, lhs -> rhs1true2false, ...]], or similar, but I thought this would have required more extensive changes to Kranc.

barrywardell commented 10 years ago

PS. The implementation is in xTensorKranc.m and could probably be ported to TensorTools relatively easily.

eschnett commented 10 years ago

The latter step can be implemented as an optimization. Common expressions can be pulled out of conditionals. They can also be pulled out of other expressions such as Max or Min, but one has to be careful there. Pulling common expressions out of sums could lead to a large performance benefit, if we can e.g. re-constitute (a+b+c)*(e+f+g) from its expanded, distributed version.