[ADOL-C] Jacobians at each time step for a PDE

[Kshitij Kulshreshtha]

Hello,

Your understanding that checkpointing will compute the final jacobian
more efficiently is quite correct. But since you need the full jacobian
in each step, there are a few possibilities.

1. you need to run the simulation again with varying number of steps
like this:

for (tstep=1; tstep == tmax; tstep++)
   trace_on(1)
   set_independants
   for (t=0; t<tstep; t++)
       model_step
   declare_dependants
   trace_off()
   jacobian(1,n_dep, n_ind, ind, J)

2. You trace only one step, but remember the previous step's full
jacobian jacobian and mutiply this steps partial jacobian with it to get
the new full jacobian

fullJ = zeromatrix
for (t=0; t<tmax; t++)
    trace_on(1)
    set_independants
    model_step
    declare_dependants
    trace_off()
    jacobian(1,n_dep, n_ind, ind, J)
    if (t == 0)
       copymatrix(fullJ, J)
    else
       fullJ = matrixmult(fullJ,J)

3. You declare every step's output as dependents on the same tape. Then
the jacobian of each step's output w.r.t input can be computed.

trace_on(1)
set_independants
for (t=0; t<tmax; t++)
    model_step
    declare_dependants
trace_off()
jacobian(1,n_dep*tmax, n_ind, ind, multiJ);


Which method is the most efficient, depends on your step implementation,
sparsity of the individual step jacobians and many other things. This is
difficult to foresee, you'll have to try it out. I would however
recommend using sparse_jac() routine instead of jacobian(), since PDE
discretisations often yield sparce matrices. As I see it, checkpointing
may not help in your case at all.


Best
Kshitij Kulshrestha

As on 2013-05-16 14:00h, Adam Poole did write:
> Hello,
> 
> I have managed to get very close to what I am trying to achieve but have
> slight problem. I have a non-linear PDE model that I am trying to
> calibrate with real data, so I wish to get a jacobian at each time step.
> If I dothe following sudo code
> 
> for (t=0; t<tmax; t++)
>     trace_on(1)
>     set_independants
>     model_step
>     declare_dependants
>     trace_off()
>     jacobian(1,n_dep, n_ind, ind, J)
> 
> Then I get jacobians at each time step, but, the tape is restarted so my
> second iteration calculates a new jacobian ignoring the infulence of the
> previous time step. The model includes several functions of the form
> x(t+1)=f(x(t),u) but this method treats x(t) as a constant.
> 
> Is there a way to get trace_on to continue the tape from the previous
> iteration if I save it? Doing the following (which is also more
> efficient with declaring inds/deps)
> 
> trace_on(1)
> set_independants
> for (t=0; t<tmax; t++)
>     model_step
> declare_dependants
> trace_off()
> jacobian(1,n_dep, n_ind, ind, J)
> 
> produces the jacobian I expect for the last time step. I had a quick
> look at the check pointing example, but am I correct in believing this
> will only produce the final jacobian in a more effcient manner. Also,
> the definition would be unweildly to implement with the way the model
> code has been written.
> 
> Many Thanks,
> Adam
> _______________________________________________
> ADOL-C mailing list
> ADOL-C at list.coin-or.org
> http://list.coin-or.org/mailman/listinfo/adol-c

-- 
Dr. Kshitij Kulshreshtha

Institut für Mathematik,
Universität Paderborn,
Warburger Straße 100,
33098 Paderborn.

Büro: A3.235

Privatanschrift:
Arnikaweg 62
33100 Paderborn.