Subject: [Paracardium-dev] ODE time integration

From: Ola Skavhaug <[email-masked]>
Date: Mon, 9 Feb 2009 10:34:33 +0100

Hi,

when solving the bidomain equations using operator splitting, we and up with
a large number of ODEs that needs to be solved at each time step.

Typically, the code will look like (excluding parameter setting and such):

ODECollection:: forward(uint n, double* v, double t, double dt)
{
   // Here n is the number of transmembrane potential dofs

   for (uint i=0; i<n; ++i) {
      states[i*ode.size()] = v[i];
      ode.forward(&states[i*ode.size()], t, dt);
      v[i] = states[i*ode.size()];
}

Hence, we have some requirements for the interface of the ode class:
1) Should be able to work with pointers (to omit allocation)
2) Should be able to integrate from time t to t+dt.
3) (Not clear from the code above) We should be able to set an initial guess
for the initial time step size

It is not all clear for me from reading the ode section of DOLFIN how to
implement these extensions.

-- 
Ola Skavhaug

---

From: Anders Logg <[email-masked]>
Date: Mon, 9 Feb 2009 12:43:40 +0100

On Mon, Feb 09, 2009 at 10:34:33AM +0100, Ola Skavhaug wrote:
> Hi,
> 
> when solving the bidomain equations using operator splitting, we and up with a
> large number of ODEs that needs to be solved at each time step.
> 

I think the ODECollection class is a good abstraction and maybe it
should even be moved to DOLFIN. Are there any reasons for not having
it in DOLFIN (like dependencies on other ODE solvers)?

> implement these extensions.

Not clear to me either. I'll think about it and add some hooks to
DOLFIN to make this work.

-- 
Anders

---

From: Ola Skavhaug <[email-masked]>
Date: Mon, 9 Feb 2009 13:12:29 +0100

> it in DOLFIN (like dependencies on other ODE solvers)?

No, if you want it in DOLFIN, I'll be happy to move it there.

> DOLFIN to make this work.

Very good. I'll think about this as well.

Ola

> --
> Anders
>
> -----BEGIN PGP SIGNATURE-----
> Version: GnuPG v1.4.9 (GNU/Linux)
>
> iEYEARECAAYFAkmQFusACgkQTuwUCDsYZdFYBACfb8zgEOY5sEXjM/j9vrPKTVY/
> dDcAnAlHddoqXcoEaqGm803OGIK8yiNG
> =/Eqx
> -----END PGP SIGNATURE-----
>
>

-- 
Ola Skavhaug

---

From: Anders Logg <[email-masked]>
Date: Mon, 9 Feb 2009 23:22:28 +0100

> Very good. I'll think about this as well.

I've pushed some stuff to allow integration from t to t + dt. Looks
like it works. Will look more at ODECollection tomorrow.

Maybe you can parallelize the assembly in the meantime and we'll be
finished by Wednesday... ;-)

-- 
Anders

---

From: Ola Skavhaug <[email-masked]>
Date: Mon, 9 Feb 2009 23:40:02 +0100

On Mon, Feb 9, 2009 at 11:22 PM, Anders Logg <[email-masked]> wrote:

> >     > when solving the bidomain equations using operator splitting, we
> and up
> >     with a
> >     > 3) (Not clear from the code above) We should be able to set an
> initial
> >     guess

That is very nice, Anders!

> finished by Wednesday... ;-)

That sounds like an excellent idea :)

I will have a look at the parallelization again tomorrow.

Wednesday we can see what interface the ODECollection abstraction should
have.

Currently, only the constructor and a forward is implemented. I'll check it
in when at work.
We need need to add support for setting parameters through the ODECollection
as well.

Ola

> iEYEARECAAYFAkmQrKQACgkQTuwUCDsYZdGqUwCcDaxWSNdoIPR4qhonOLbtBj61
> PyMAmwVYdzpcW20aoqNOBaDuRl1v1Qsb
> =3j8Q

-- 
Ola Skavhaug

---

From: Anders Logg <[email-masked]>
Date: Tue, 10 Feb 2009 08:10:28 +0100

On Mon, Feb 09, 2009 at 11:40:02PM +0100, Ola Skavhaug wrote:
> Wednesday we can see what interface the ODECollection abstraction should have.
> 
> Currently, only the constructor and a forward is implemented. I'll check it in
> when at work.
> We need need to add support for setting parameters through the ODECollection as
> well.

I'm about to implement the ODECollection class in DOLFIN. I think I
have a good idea how to implement it but you might want to change it
after I've made the first iteration.

-- 
Anders

---

From: Ola Skavhaug <[email-masked]>
Date: Tue, 10 Feb 2009 08:48:48 +0100

On Tue, Feb 10, 2009 at 8:10 AM, Anders Logg <[email-masked]> wrote:

> >     >     > Typically, the code will look like (excluding parameter
> setting and
> >     >     such):
> >     >     >
> >     >     > It is not all clear for me from reading the ode section of
> DOLFIN
> >     how to

Fine. Lets hope we're more or less in sync when it comes to this
abstraction.

Ola

> iEYEARECAAYFAkmRKGQACgkQTuwUCDsYZdE7FgCdEMRitfnVlG7aiYlVoRv+my6W
> nsYAnjJmJMqFkR1UYsmmJSKMHJrbXKVl
> =E2wa

-- 
Ola Skavhaug

---

From: Anders Logg <[email-masked]>
Date: Tue, 10 Feb 2009 14:48:45 +0100

On Tue, Feb 10, 2009 at 08:48:48AM +0100, Ola Skavhaug wrote:
> Fine. Lets hope we're more or less in sync when it comes to this abstraction.
>  
> Ola

Let's hope so. I've pushed it to DOLFIN. See what you think. There's a
demo in

  demo/ode/collection/cpp/

-- 
Anders

---

From: Ola Skavhaug <[email-masked]>
Date: Tue, 10 Feb 2009 15:13:52 +0100

On Tue, Feb 10, 2009 at 2:48 PM, Anders Logg <[email-masked]> wrote:

> >     >     > On Mon, Feb 9, 2009 at 12:43 PM, Anders Logg <[email-masked]
> >
> >     wrote:
> >     >     >
> >     >     >     > when solving the bidomain equations using operator
> splitting,
> >     we
> >     >     and up
> >     >     >     with a
> >     >     >     > large number of ODEs that needs to be solved at each
> time
> >     step.
> >     >     >     >
> >     >     >     > ODECollection:: forward(uint n, double* v, double t,
> double
> >     dt)
> >     >     >
> >     >     >     I think the ODECollection class is a good abstraction and
> maybe
> >     it
> >     >     >     should even be moved to DOLFIN. Are there any reasons for
> not
> >     having
> >     >     >     > Hence, we have some requirements for the interface of
> the ode
> >     >     class:
> >     >     >     > 3) (Not clear from the code above) We should be able to
> set
> >     an
> >     >     initial
> >     >     >     guess
> >     >     >     Not clear to me either. I'll think about it and add some
> hooks
> >     to
> >     > Currently, only the constructor and a forward is implemented. I'll
> check
> >     it in
> >     > when at work.

This looks very nice. I have to run some profiling before I can say anything
about the efficiency, though.

The next point to consider is the possible need to vary certain parameters
from one ODE to another.

In PyCC, we have implemented an efficient way of getting and setting
parameters for ODEs. I guess a similar construction is possible for ODEs,
being children of Parametrized.

Baiscally, we have a few methods for doing this:

* int numParameters() const;
* void setParameter(int idx, double value);
* void setParameters(const double* values);
* double getParameter(int idx) const;
* void getParameters(double* values) const;
* std::string getParameterName(int idx) const;

The nice thing about this construction is that we can send a parameter
vector to ODECollection and use this vector to set various parameters during
the ode solve loop. One important parameter that is typically set during
this loop is the stimulus current.

Ola

> iEYEARECAAYFAkmRhbwACgkQTuwUCDsYZdH3wQCfYI1hs9Y4x7BjPT6V3x6h5HI1
> iF0An1F/cOJ5MIRnBn9iSDcpLrpUi7K2
> =wyg8

-- 
Ola Skavhaug

---

From: Anders Logg <[email-masked]>
Date: Tue, 10 Feb 2009 19:23:32 +0100

> about the efficiency, though.

Yes, I think a lot can be done but we need profiling to know where to look.

> The next point to consider is the possible need to vary certain parameters from
> one ODE to another.
> 
> In PyCC, we have implemented an efficient way of getting and setting parameters
> for ODEs. I guess a similar construction is possible for ODEs, being children
> of Parametrized.
> 
> The nice thing about this construction is that we can send a parameter vector
> to ODECollection and use this vector to set various parameters during the ode
> solve loop. One important parameter that is typically set during this loop is
> the stimulus current.

Sounds like it would be better to have a callback so a user is allowed
to do anything necessary with the ODE in between calls.

I've added a callback to ODECollection that may optionally be
overloaded to change parameters:

  /// Optional user-defined update, called between solves
  virtual void update(ODE& ode, real* u, real t, uint system);

-- 
Anders