matt-weinstein / adigator

Matlab Algorithmic Differentiation Toolbox
25 stars 5 forks source link

Variable reuse in user function leads to non-executable Jacobian files #13

Open rschiemann1 opened 4 years ago

rschiemann1 commented 4 years ago

I have a right hand side function according to the following definition:

y' = f(t, y, p)
where: 
   y' is the time derivative of the solution vector y
   t is the time
   y is the solution vector
   p is a structure holding numeric parameters.

In my code, the function f is called 'sysfun', the parameter structure p is called 'Params'.

As is usual in Matlab ODE solvers, the system function to be integrated needs to give the vector of first order time derivatives as output (compare to functional form at the beginning of this post). In my code, I do the following in numerous places:

kk = 1; 
dy(IndexMap1(kk)) = some code involving y(kk) and y(kk+1).
kk = n; 
dy(IndexMap1(kk)) = some other code involving y(kk) and y(kk-1).
kk = 2:n-1
dy(IndexMap1(kk)) = some code involving y(kk) and y(kk+1) or y(kk-1).

As my be apparent, such code results from spatial discretization of partial differential integrations into n discrete volumes and then using the Matlab ODE solvers to integrate the resulting set of ODEs over time.

The variable kk is reused for a few different discretized PDEs, which seems to be the root of the problem I am about to report.

The issue that I now have is the following:

Creating the Jacobian file with the following lines of code finishes successfully without any errors or warnings:

ay = adigatorCreateDerivInput(size(y0),'y');
adiOpts = adigatorOptions();
adiOpts.overwrite = 1;
output = adigatorGenJacFile('sysfun',{1,ay,Params},adiOpts);

In above code, y0 is a vector of initial values of the system.

However, when I then try to run the generated Jacobian file (e.g. with the command sysfun_Jac(1,y0,Params)), I get the following error message:

Unable to perform assignment because dot indexing is not supported for variables of this type.

The reason is that the generated Jacobian code also reuses the variable kk, which itself is not a problem. However, for example the user line kk = 1 is transformed to kk = 1 in one place and to kk.f = 1 in some other please. I do not see a reason why adigator decides to do so (the user codes of the two cases are quite similar to my eye), but using kk as a scalar variable first and then assigning a field to it as if it were a structure is of course not possible.

I was able to work around this issue by not reusing kk in my user code over and over, but use different variables names, i.e. kk1, kk2, kk3, ….