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

Allow selection of finite differencing scheme for a calculation at runtime #1

Open ianhinder opened 14 years ago

ianhinder commented 14 years ago

It is currently not possible to specify the finite differencing operators used by a calculation at runtime. At the moment, one can use Mathematica to generate multiple copies of a calculation, each with a different set of FD operators (i.e. different FD schemes). This generates multiple source files which all need to be compiled. This is inelegant and takes a lot of time during compilation.

We could instead have a conditional at the start of the calculation which selects which FD operators to use to compute each derivative. This was not done before because we thought there might be a performance hit from a conditional at every point in the loop, but this might not be an issue on modern systems, especially if you can indicate to the compiler that this variable does not change in some sense.

In a calculation, we currently have "named derivatives". A finite differencing "scheme" would be a mapping of these named derivatives to finite differencing operators. We could give each scheme a name, and the user specifies the mapping. We could then allow the user to specify a scheme for the whole thorn in the parameter file. Additionally, this could be overridden on a per-calculation basis, for example to use a different scheme for some calculations.

ianhinder commented 12 years ago

Instead of having named schemes, I have made it possible for finite difference operators to depend on integer parameters (see 9737d7a9f4e526f4aa14d9701c255ab06789f445). This addresses the most immediate need, which is to allow the finite difference order to be selectable at run-time.

One idea for implementing different differencing schemes would be to have an option

PartialDerivativeSchemes -> {{name, defs}, {name, defs}, ...}

where defs is what we currently provide in PartialDerivatives. name could be "centered", "upwinded", etc.

ianhinder commented 12 years ago

The complication of having named schemes is that we would need to do a string comparison. The result of this would have to be computed outside the grid function loop and mapped to an integer comparison. You would then need a set of multi-level switch statements in the inner loop.

Another way to approach it would be to allow difference operators to depend on more than one parameter, including keyword parameters. This would be easier to incorporate in the current code. So you could have something like

PartialDerivatives -> {
  PDAdvect[i_] ->
    Switch[advectScheme,
      "upwind", UpwindDerivative[beta,fdOrder,i],
      "upwind-split", Sign[beta[i]] UpwindAntisymmetricDerivative[fdOrder,i] + 
        Abs[beta[i]]UpwindSymmetricDerivative[fdOrder,i],
      "centered", CenteredDerivative[fdOrder,i]]
}

This might be preferable to having multiple sets of PartialDerivative definitions, as it means that different schemes don't need to duplicate operators that they have in common.