1D PDE solver for Octave

classic Classic list List threaded Threaded
24 messages Options
12
Reply | Threaded
Open this post in threaded view
|

1D PDE solver for Octave

Bill Greene-3
I'm working on a 1D PDE solver for Octave that is somewhat
similar to the MATLAB function pdepe.

I would like to make this available to the general Octave
community but am not sure what is the best strategy for this.
So, I'm soliciting suggestions.

The function is implemented as a C++ mex file. It depends very
heavily on the Sundials ida ODE solver, somewhat less heavily
on the Eigen C++ matrix class library, and a couple of libraries
from Boost.

So, my main issue is how to deal with these three dependencies 
when distributing the code? I suppose I could build binary
versions myself for a few platforms and package these with the
source. I haven't come across anyone else distributing Octave
functions this way. I perused a few projects on octave-forge
but haven't seen any that have faced this same issue; the pkg
install capability is certainly convenient for users.

I know that some work is ongoing using Sundials ida to improve the
Octave ODE solvers. But I'm unclear on when or if this might result
in the Sundials libraries being available in Octave.

So, as I say, I'll appreciate any suggestions on how to proceed.

Bill Greene
Reply | Threaded
Open this post in threaded view
|

Re: 1D PDE solver for Octave

Carlo de Falco-2

On 29 Jun 2016, at 19:29, Bill Greene <[hidden email]> wrote:

> The function is implemented as a C++ mex file. It depends very
> heavily on the Sundials ida ODE solver, somewhat less heavily
> on the Eigen C++ matrix class library, and a couple of libraries
> from Boost.
>

Hi,

Thanks for your interest in contributing to Octave!

I'd like to know more about what you are doing in order
to provide more informed suggestions.

I'm not sure what you intend by "C++ MEX file" MEX is a C API,
Octave has its own C++ API.

You say you are relying on sundias IDA, are you using
IDA's own C API or the MEX interface?

There is another ongoing project to implement matlab
compatible ode15s and ode15i via and interface to sundias IDA [1,2].

So in the near future there will be also an oct-file interface to
IDA which based on initial benchamarking is substantially
faster than the MEX file interface distributed with sundials.

I am also curious about what algorithm does your code exactly
implement because IDA and PDEPE are based on completely different
mathematical methods so I'm not sure what you are actually doing ...

c.


[1] http://gsoc2016ode15s.blogspot.it/
[2] https://bitbucket.org/Francesco_Faccio/octave/

Reply | Threaded
Open this post in threaded view
|

Re: 1D PDE solver for Octave

Bill Greene-3
Thanks for taking the time to reply although I am disappointed
that you didn't address any of my issues.

I am using the C API to ida; my code is written in C++.

My code is currently using a simple linear interpolation in the
spatial dimension and then using ida to solve the resulting system
of ODE. I intend to improve the spatial discretization in the future
but the current implementation could certainly be considered "inferior"
to pdepe. I'm not trying to duplicate pdepe-- just support the same
problem class and the command line interface as closely as possible.

I should say one more thing about Sundials ida. The latest version of 
ida is only just "adequate" to support my current implementation,
mainly because its sparse matrix support is rather rudimentary. I am
hoping they will improve ida in this area and expect that for the
foreseeable future I'll want to be using the very latest ida release.

I am aware of the ongoing Octave ode15s work; I alluded to that in my original
message.

I would be happy for anyone to try the current code. At this point, I
consider it to be solid beta-quality code.

Thanks.

Bill

On Thu, Jun 30, 2016 at 8:31 AM, Carlo De Falco <[hidden email]> wrote:

On 29 Jun 2016, at 19:29, Bill Greene <[hidden email]> wrote:

> The function is implemented as a C++ mex file. It depends very
> heavily on the Sundials ida ODE solver, somewhat less heavily
> on the Eigen C++ matrix class library, and a couple of libraries
> from Boost.
>

Hi,

Thanks for your interest in contributing to Octave!

I'd like to know more about what you are doing in order
to provide more informed suggestions.

I'm not sure what you intend by "C++ MEX file" MEX is a C API,
Octave has its own C++ API.

You say you are relying on sundias IDA, are you using
IDA's own C API or the MEX interface?

There is another ongoing project to implement matlab
compatible ode15s and ode15i via and interface to sundias IDA [1,2].

So in the near future there will be also an oct-file interface to
IDA which based on initial benchamarking is substantially
faster than the MEX file interface distributed with sundials.

I am also curious about what algorithm does your code exactly
implement because IDA and PDEPE are based on completely different
mathematical methods so I'm not sure what you are actually doing ...

c.


[1] http://gsoc2016ode15s.blogspot.it/
[2] https://bitbucket.org/Francesco_Faccio/octave/

Reply | Threaded
Open this post in threaded view
|

Re: 1D PDE solver for Octave

Carlo de Falco-2

On 30 Jun 2016, at 15:08, Bill Greene <[hidden email]> wrote:

> Thanks for taking the time to reply although I am disappointed
> that you didn't address any of my issues.

As I said in my previus email I needed more info in order to address
some of them.

> I am using the C API to ida; my code is written in C++.

This is the same approach Franceso is using in ode15i/ode15s.

There may be lots of code duplication, it would be nice to compare
your implementations ...

> My code is currently using a simple linear interpolation in the
> spatial dimension and then using ida to solve the resulting system
> of ODE.

So to summarize, given a parabolic PDE you discretize in space
using linear finite elements / first order finite differences
the solve the resulting ODE with BDF time-stepping via IDA,
is this correct?

> I intend to improve the spatial discretization in the future
> but the current implementation could certainly be considered "inferior"
> to pdepe. I'm not trying to duplicate pdepe-- just support the same
> problem class and the command line interface as closely as possible.

If my interpretation of your algorithm above is correct, it sounds
like, in order to reduce code duplication, this could be implemented
like a thin wrapper around BIM and ode15i or daspk, but we need to see
the code to be more specific.

> I should say one more thing about Sundials ida. The latest version of
> ida is only just "adequate" to support my current implementation,
> mainly because its sparse matrix support is rather rudimentary. I am
> hoping they will improve ida in this area and expect that for the
> foreseeable future I'll want to be using the very latest ida release.
>
> I am aware of the ongoing Octave ode15s work; I alluded to that in my original
> message.

I missed that, sorry. Replying below.

>> I know that some work is ongoing using Sundials ida to improve the
>> Octave ODE solvers. But I'm unclear on when or if this might result
>> in the Sundials libraries being available in Octave.

The changes for linking Octave to Sundials in Francesco's code are essentially
done, but I'm not sure whether they'll be merged into core for the 4.2 or 4.4
release ...

If you want to start using those you can just clone Francesco's repo
from bitbucket and build Octave from there.

>
> I would be happy for anyone to try the current code. At this point, I
> consider it to be solid beta-quality code.

sure, were is the code?

>
> Thanks.
>
> Bill


c.
Reply | Threaded
Open this post in threaded view
|

Re: 1D PDE solver for Octave

Bill Greene-3
>sure, were is the code?
OK, great.

The source code is here:

Unfortunately, my approach to building this is somewhat ad-hoc
at this point so you might find that challenging. I can send you a
Linux Makefile if you are sufficiently motivated.

Alternatively, I have a pre-built Linux binary here:


that seems to work fine with Octave 4.02 on my system.

Let me know if you run into problems with either approach.

Thanks.

Bill

On Thu, Jun 30, 2016 at 9:54 AM, Carlo De Falco <[hidden email]> wrote:

On 30 Jun 2016, at 15:08, Bill Greene <[hidden email]> wrote:

> Thanks for taking the time to reply although I am disappointed
> that you didn't address any of my issues.

As I said in my previus email I needed more info in order to address
some of them.

> I am using the C API to ida; my code is written in C++.

This is the same approach Franceso is using in ode15i/ode15s.

There may be lots of code duplication, it would be nice to compare
your implementations ...

> My code is currently using a simple linear interpolation in the
> spatial dimension and then using ida to solve the resulting system
> of ODE.

So to summarize, given a parabolic PDE you discretize in space
using linear finite elements / first order finite differences
the solve the resulting ODE with BDF time-stepping via IDA,
is this correct?

> I intend to improve the spatial discretization in the future
> but the current implementation could certainly be considered "inferior"
> to pdepe. I'm not trying to duplicate pdepe-- just support the same
> problem class and the command line interface as closely as possible.

If my interpretation of your algorithm above is correct, it sounds
like, in order to reduce code duplication, this could be implemented
like a thin wrapper around BIM and ode15i or daspk, but we need to see
the code to be more specific.

> I should say one more thing about Sundials ida. The latest version of
> ida is only just "adequate" to support my current implementation,
> mainly because its sparse matrix support is rather rudimentary. I am
> hoping they will improve ida in this area and expect that for the
> foreseeable future I'll want to be using the very latest ida release.
>
> I am aware of the ongoing Octave ode15s work; I alluded to that in my original
> message.

I missed that, sorry. Replying below.

>> I know that some work is ongoing using Sundials ida to improve the
>> Octave ODE solvers. But I'm unclear on when or if this might result
>> in the Sundials libraries being available in Octave.

The changes for linking Octave to Sundials in Francesco's code are essentially
done, but I'm not sure whether they'll be merged into core for the 4.2 or 4.4
release ...

If you want to start using those you can just clone Francesco's repo
from bitbucket and build Octave from there.

>
> I would be happy for anyone to try the current code. At this point, I
> consider it to be solid beta-quality code.

sure, were is the code?

>
> Thanks.
>
> Bill


c.

Reply | Threaded
Open this post in threaded view
|

Re: 1D PDE solver for Octave

Juan Pablo Carbajal-2
On Thu, Jun 30, 2016 at 7:56 PM, Bill Greene <[hidden email]> wrote:

>>sure, were is the code?
> OK, great.
>
> The source code is here:
> https://github.com/wgreene310/pde1d
>
> Unfortunately, my approach to building this is somewhat ad-hoc
> at this point so you might find that challenging. I can send you a
> Linux Makefile if you are sufficiently motivated.
>
You should distribute your code with a makefile if you have the
intention that other people test it.


> Alternatively, I have a pre-built Linux binary here:
>
> https://drive.google.com/folderview?id=0B0HVh_MUHujMU0MzY2FxSlVUWWs&usp=sharing
>
> that seems to work fine with Octave 4.02 on my system.
>
> Let me know if you run into problems with either approach.
>
> Thanks.
>
> Bill
>
> On Thu, Jun 30, 2016 at 9:54 AM, Carlo De Falco <[hidden email]>
> wrote:
>>
>>
>> On 30 Jun 2016, at 15:08, Bill Greene <[hidden email]> wrote:
>>
>> > Thanks for taking the time to reply although I am disappointed
>> > that you didn't address any of my issues.
>>
>> As I said in my previus email I needed more info in order to address
>> some of them.
>>
>> > I am using the C API to ida; my code is written in C++.
>>
>> This is the same approach Franceso is using in ode15i/ode15s.
>>
>> There may be lots of code duplication, it would be nice to compare
>> your implementations ...
>>
>> > My code is currently using a simple linear interpolation in the
>> > spatial dimension and then using ida to solve the resulting system
>> > of ODE.
>>
>> So to summarize, given a parabolic PDE you discretize in space
>> using linear finite elements / first order finite differences
>> the solve the resulting ODE with BDF time-stepping via IDA,
>> is this correct?
>>
>> > I intend to improve the spatial discretization in the future
>> > but the current implementation could certainly be considered "inferior"
>> > to pdepe. I'm not trying to duplicate pdepe-- just support the same
>> > problem class and the command line interface as closely as possible.
>>
>> If my interpretation of your algorithm above is correct, it sounds
>> like, in order to reduce code duplication, this could be implemented
>> like a thin wrapper around BIM and ode15i or daspk, but we need to see
>> the code to be more specific.
>>
>> > I should say one more thing about Sundials ida. The latest version of
>> > ida is only just "adequate" to support my current implementation,
>> > mainly because its sparse matrix support is rather rudimentary. I am
>> > hoping they will improve ida in this area and expect that for the
>> > foreseeable future I'll want to be using the very latest ida release.
>> >
>> > I am aware of the ongoing Octave ode15s work; I alluded to that in my
>> > original
>> > message.
>>
>> I missed that, sorry. Replying below.
>>
>> >> I know that some work is ongoing using Sundials ida to improve the
>> >> Octave ODE solvers. But I'm unclear on when or if this might result
>> >> in the Sundials libraries being available in Octave.
>>
>> The changes for linking Octave to Sundials in Francesco's code are
>> essentially
>> done, but I'm not sure whether they'll be merged into core for the 4.2 or
>> 4.4
>> release ...
>>
>> If you want to start using those you can just clone Francesco's repo
>> from bitbucket and build Octave from there.
>>
>> >
>> > I would be happy for anyone to try the current code. At this point, I
>> > consider it to be solid beta-quality code.
>>
>> sure, were is the code?
>>
>> >
>> > Thanks.
>> >
>> > Bill
>>
>>
>> c.
>
>

Reply | Threaded
Open this post in threaded view
|

Re: 1D PDE solver for Octave

Bill Greene-3
I just checked in a Makefile that works for me on linux.

But my strong recommendation, at this point, is just to use the pre-built
shared library I provided the link to.

As you recall, the main point of my original post is how to deal with these build
issues in a reasonable way. My belief is that very few Octave users will try this
pde solver if they have to build it themselves.


On Mon, Jul 4, 2016 at 4:30 AM, Juan Pablo Carbajal <[hidden email]> wrote:
On Thu, Jun 30, 2016 at 7:56 PM, Bill Greene <[hidden email]> wrote:
>>sure, were is the code?
> OK, great.
>
> The source code is here:
> https://github.com/wgreene310/pde1d
>
> Unfortunately, my approach to building this is somewhat ad-hoc
> at this point so you might find that challenging. I can send you a
> Linux Makefile if you are sufficiently motivated.
>
You should distribute your code with a makefile if you have the
intention that other people test it.


> Alternatively, I have a pre-built Linux binary here:
>
> https://drive.google.com/folderview?id=0B0HVh_MUHujMU0MzY2FxSlVUWWs&usp=sharing
>
> that seems to work fine with Octave 4.02 on my system.
>
> Let me know if you run into problems with either approach.
>
> Thanks.
>
> Bill
>
> On Thu, Jun 30, 2016 at 9:54 AM, Carlo De Falco <[hidden email]>
> wrote:
>>
>>
>> On 30 Jun 2016, at 15:08, Bill Greene <[hidden email]> wrote:
>>
>> > Thanks for taking the time to reply although I am disappointed
>> > that you didn't address any of my issues.
>>
>> As I said in my previus email I needed more info in order to address
>> some of them.
>>
>> > I am using the C API to ida; my code is written in C++.
>>
>> This is the same approach Franceso is using in ode15i/ode15s.
>>
>> There may be lots of code duplication, it would be nice to compare
>> your implementations ...
>>
>> > My code is currently using a simple linear interpolation in the
>> > spatial dimension and then using ida to solve the resulting system
>> > of ODE.
>>
>> So to summarize, given a parabolic PDE you discretize in space
>> using linear finite elements / first order finite differences
>> the solve the resulting ODE with BDF time-stepping via IDA,
>> is this correct?
>>
>> > I intend to improve the spatial discretization in the future
>> > but the current implementation could certainly be considered "inferior"
>> > to pdepe. I'm not trying to duplicate pdepe-- just support the same
>> > problem class and the command line interface as closely as possible.
>>
>> If my interpretation of your algorithm above is correct, it sounds
>> like, in order to reduce code duplication, this could be implemented
>> like a thin wrapper around BIM and ode15i or daspk, but we need to see
>> the code to be more specific.
>>
>> > I should say one more thing about Sundials ida. The latest version of
>> > ida is only just "adequate" to support my current implementation,
>> > mainly because its sparse matrix support is rather rudimentary. I am
>> > hoping they will improve ida in this area and expect that for the
>> > foreseeable future I'll want to be using the very latest ida release.
>> >
>> > I am aware of the ongoing Octave ode15s work; I alluded to that in my
>> > original
>> > message.
>>
>> I missed that, sorry. Replying below.
>>
>> >> I know that some work is ongoing using Sundials ida to improve the
>> >> Octave ODE solvers. But I'm unclear on when or if this might result
>> >> in the Sundials libraries being available in Octave.
>>
>> The changes for linking Octave to Sundials in Francesco's code are
>> essentially
>> done, but I'm not sure whether they'll be merged into core for the 4.2 or
>> 4.4
>> release ...
>>
>> If you want to start using those you can just clone Francesco's repo
>> from bitbucket and build Octave from there.
>>
>> >
>> > I would be happy for anyone to try the current code. At this point, I
>> > consider it to be solid beta-quality code.
>>
>> sure, were is the code?
>>
>> >
>> > Thanks.
>> >
>> > Bill
>>
>>
>> c.
>
>

Reply | Threaded
Open this post in threaded view
|

Re: 1D PDE solver for Octave

Carlo de Falco-2
In reply to this post by Carlo de Falco-2

On 30 Jun 2016, at 15:54, Carlo De Falco <[hidden email]> wrote:

>> I am using the C API to ida; my code is written in C++.
>
> This is the same approach Franceso is using in ode15i/ode15s.


I stand corrected.

Looking at your code I see it is in C++ but it is using the Matlab compatible
MEX API (which is a C API).

This is not what Francesco has been doing.
Francesco is using the Octave C++ API:

https://www.gnu.org/software/octave/doc/interpreter/Oct_002dFiles.html#Oct_002dFiles

"Oct-files are pieces of C++ code that have been compiled with the Octave API into a dynamically loadable object. They take their name from the file which contains the object which has the extension .oct."

The latter is known to be much more efficient.

Apart from that there is not much more I can comment.

Actually whithout a working Makefile and whithout one single line of documentation,
not even a docstring describing input syntax, I sincerely doubt you will ever find anyone
willing to spend more time looking at your code.

Can you provide a least one usage example?
What class of problems can your code solve?
How would you invoke your function?


For example, the attached file shows how to add a docstring and a demo
to your function. Type "help parabolic1d" to see the help and type
"demo parabolic1d" to see the code of the demo and the resulting plot.


c.



Reply | Threaded
Open this post in threaded view
|

Re: 1D PDE solver for Octave

Jordi Gutiérrez Hermoso-2
In reply to this post by Bill Greene-3
On Mon, 2016-07-04 at 08:57 -0400, Bill Greene wrote:
> But my strong recommendation, at this point, is just to use the pre-built
> shared library I provided the link to.

It appears that you are unfamiliar with cross-platform binary
compatibility. The mex file you built is linked against specific
versions of libraries that only exist in very specific systems. It's
also built for 32-bit Intel Linux, which will not work for other
architectures. Many of us thus can't use your codex without building it
ourselves.

> My belief is that very few Octave users will try this pde solver if
> they have to build it themselves.

You're currently talking to people who are used to building Octave.
Later we can worry about packaging and distributing your built code or
we can let our downstream packagers do it for us.

Thanks for providing a Makefile. It's a starting point.

- Jordi G. H.







Reply | Threaded
Open this post in threaded view
|

Re: 1D PDE solver for Octave

Carlo de Falco-2
In reply to this post by Carlo de Falco-2

On 5 Jul 2016, at 00:58, Carlo De Falco <[hidden email]> wrote:

>
> For example, the attached file shows how to add a docstring and a demo
> to your function. Type "help parabolic1d" to see the help and type
> "demo parabolic1d" to see the code of the demo and the resulting plot.

oops, the attachment was missing!
c.


parabolic1d.m (2K) Download Attachment
Reply | Threaded
Open this post in threaded view
|

Re: 1D PDE solver for Octave

Bill Greene-3
In reply to this post by Carlo de Falco-2
I had included a help file, pde1d.m, in the binary distribution.
I just uploaded this file to https://github.com/wgreene310/pde1d
I included a "demo" at the bottom of this file hoping that the demo
function would find it there; but "demo pde1d" does not work for me.
I looked around for documentation or an example of how to create
a demo for a so-called "builtin" function but could not find one-- e.g.
I see that "demo lsode" doesn't do anything useful.

If there is a trick to doing this, please let me know.

>What class of problems can your code solve?
>How would you invoke your function?
Probably I was a bit vague about this function relative to the MATLAB
pdepe function. As far as problem class, the two main limitations that
I know of compared with pdepe are that it doesn't support complex
coefficients (yet) and does not support the "events" option of pdepe.
Beyond those two limitations, if anyone has an example that "works"
in pdepe but fails in pde1d, I would very much like to see it.

I am currently working on better documentation but it is not yet to the
first draft stage. 

On Mon, Jul 4, 2016 at 6:58 PM, Carlo De Falco <[hidden email]> wrote:

On 30 Jun 2016, at 15:54, Carlo De Falco <[hidden email]> wrote:

>> I am using the C API to ida; my code is written in C++.
>
> This is the same approach Franceso is using in ode15i/ode15s.


I stand corrected.

Looking at your code I see it is in C++ but it is using the Matlab compatible
MEX API (which is a C API).

This is not what Francesco has been doing.
Francesco is using the Octave C++ API:

https://www.gnu.org/software/octave/doc/interpreter/Oct_002dFiles.html#Oct_002dFiles

"Oct-files are pieces of C++ code that have been compiled with the Octave API into a dynamically loadable object. They take their name from the file which contains the object which has the extension .oct."

The latter is known to be much more efficient.

Apart from that there is not much more I can comment.

Actually whithout a working Makefile and whithout one single line of documentation,
not even a docstring describing input syntax, I sincerely doubt you will ever find anyone
willing to spend more time looking at your code.

Can you provide a least one usage example?
What class of problems can your code solve?
How would you invoke your function?


For example, the attached file shows how to add a docstring and a demo
to your function. Type "help parabolic1d" to see the help and type
"demo parabolic1d" to see the code of the demo and the resulting plot.


c.



Reply | Threaded
Open this post in threaded view
|

Re: 1D PDE solver for Octave

Bill Greene-3
In reply to this post by Jordi Gutiérrez Hermoso-2
>It's also built for 32-bit Intel Linux, which will not work for other architectures
Of course.
I built the mex file for Octave 4.0 and lnux32 assuming that:
a) it would run on any 32-bit linux that Octave 4.0 runs on 
b) most Octave developers/maintainers are running on this platform

If there are things I can do to build the .so so it runs on more linux32 systems,
I would like to hear them (this is really at the heart of my original request for
help).

On Mon, Jul 4, 2016 at 8:04 PM, Jordi Gutiérrez Hermoso <[hidden email]> wrote:
On Mon, 2016-07-04 at 08:57 -0400, Bill Greene wrote:
> But my strong recommendation, at this point, is just to use the pre-built
> shared library I provided the link to.

It appears that you are unfamiliar with cross-platform binary
compatibility. The mex file you built is linked against specific
versions of libraries that only exist in very specific systems. It's
also built for 32-bit Intel Linux, which will not work for other
architectures. Many of us thus can't use your codex without building it
ourselves.

> My belief is that very few Octave users will try this pde solver if
> they have to build it themselves.

You're currently talking to people who are used to building Octave.
Later we can worry about packaging and distributing your built code or
we can let our downstream packagers do it for us.

Thanks for providing a Makefile. It's a starting point.

- Jordi G. H.







Reply | Threaded
Open this post in threaded view
|

Re: 1D PDE solver for Octave

Carlo de Falco-2
In reply to this post by Bill Greene-3

On 5 Jul 2016, at 15:34, Bill Greene <[hidden email]> wrote:

> but "demo pde1d" does not work for me.

I haven't tried it for pde1d.m because I was not able to build pde1d yet
but it works for any mex file I tried:

- let function "fun" be defined in file "fun.mex"
- let file "fun.m" contain a docstring and a demo block
- let both be in the path

typing "help fun" displays the docstring from "fun.m"
typing "demo fun" runs the demo from "fun.m"

> I looked around for documentation or an example of how to create
> a demo for a so-called "builtin" function but could not find one-- e.g.

pde1d is not a builtin function, it is a MEX file, so I don't think
this is relevant for the case at hand.

Anyway, if you have a builtin function defined in a C++ file named
"pippo.cc" and "pippo.cc" is in your path, typing

 demo pippo

or

 demo pippo.cc

will  run the demo blocks inside "pippo.cc"

You can try this with the attached file.

> I see that "demo lsode" doesn't do anything useful.

that's because maybe lsode has no demos?

c.


pippo.cc (54 bytes) Download Attachment
Reply | Threaded
Open this post in threaded view
|

Re: 1D PDE solver for Octave

Carlo de Falco-2
In reply to this post by Bill Greene-3

On 5 Jul 2016, at 15:42, Bill Greene <[hidden email]> wrote:

> b) most Octave developers/maintainers are running on this platform
Actually, I believe 64bit systems are more common nowadays, aren't they?
c.

Reply | Threaded
Open this post in threaded view
|

Re: 1D PDE solver for Octave

Carlo de Falco-2
In reply to this post by Bill Greene-3

On 4 Jul 2016, at 14:57, Bill Greene <[hidden email]> wrote:

> I just checked in a Makefile that works for me on linux.

I rearranged your code and fixed your makefile so that it works
for me on OSX too.

This should work with minimal editing on any other platform with
Eigen, Octave and Sundials installed too.

Your demo does not work because its syntax is wrong.

Only in the first line you should have "%!demo" without
a space after "%!".

Even if I fix this the demo does nothing useful, it defines a
bunch of functions but doesn't run anything.

c.


pde1d-master.tgz (19K) Download Attachment
Reply | Threaded
Open this post in threaded view
|

Re: 1D PDE solver for Octave

Bill Greene-3
Glad to hear you were able to build on OSX. I don't currently have access to that platform.

Regarding the demo, it looks to me like the first "demo" line is
%!demo
I don't see an extra space. I haven't yet tried on linux but when I do
demo pde1d
on windows, I get this error:
error: element number 1 undefined in return list
error: called from
    demo at line 111 column 13

>Even if I fix this the demo does nothing useful, it defines a
>bunch of functions but doesn't run anything.

OK, I just took the contents of one of my test m-files and put
a %! in front of every line. So maybe there should be
%!heatCond
as the last line in the demo section to call the "main" function?


On Tue, Jul 5, 2016 at 11:42 AM, Carlo De Falco <[hidden email]> wrote:

On 4 Jul 2016, at 14:57, Bill Greene <[hidden email]> wrote:

> I just checked in a Makefile that works for me on linux.

I rearranged your code and fixed your makefile so that it works
for me on OSX too.

This should work with minimal editing on any other platform with
Eigen, Octave and Sundials installed too.

Your demo does not work because its syntax is wrong.

Only in the first line you should have "%!demo" without
a space after "%!".

Even if I fix this the demo does nothing useful, it defines a
bunch of functions but doesn't run anything.

c.


Reply | Threaded
Open this post in threaded view
|

Re: 1D PDE solver for Octave

Carlo de Falco-2
Bill,

On 5 Jul 2016, at 18:31, Bill Greene <[hidden email]> wrote:

> Glad to hear you were able to build on OSX. I don't currently have access to that platform.

This version of the Makefile should be relatively platform independent
as it relies on the "mkoctfile".

This is the standard approach used in all Octave packages.
Using a slightly improved version of this file you should
be able to quite easily reorganize your code into a package
using thie procedure described here:

https://www.gnu.org/software/octave/doc/interpreter/Creating-Packages.html#Creating-Packages

this is how Octave extesions are

> Regarding the demo, it looks to me like the first "demo" line is
> %!demo
> I don't see an extra space.

there shoudn't be any, the first line must be exactly

%!demo

the space should be in the other lines, as described here:

https://www.gnu.org/software/octave/doc/interpreter/Demonstration-Functions.html#Demonstration-Functions


> I haven't yet tried on linux but when I do
> demo pde1d
> on windows, I get this error:
> error: element number 1 undefined in return list
> error: called from
>     demo at line 111 column 13
>
> >Even if I fix this the demo does nothing useful, it defines a
> >bunch of functions but doesn't run anything.
>
> OK, I just took the contents of one of my test m-files and put
> a %! in front of every line. So maybe there should be
> %!heatCond
> as the last line in the demo section to call the "main" function?

typing "demo pde1d" is equivalent to copying the contents of the demo
blocks and typing "enter". If you do that, do you get any interesting result?

In the manual page for demo linked above, please read carefully the paragraph about
defining functions in a demo block:

'... Finally, because demo evaluates within a function context it is not possible to define new functions within the code. Anonymous functions make a good substitute in most instances. If function blocks must be used then the code eval (example ("function", n)) will allow Octave to see them. This has its own problems, however, as eval only evaluates one line or statement at a time. In this case the function declaration must be wrapped with "if 1 <demo stuff> endif" where "if" is on the same line as "demo" ...'

c.



Reply | Threaded
Open this post in threaded view
|

Re: 1D PDE solver for Octave

Bill Greene-3
>This version of the Makefile should be relatively platform independent
>as it relies on the "mkoctfile".
Yes, I see the value of that.
Nevertheless, I had to make quite a few changes to your Makefile in order
to build on debian linux. The pre-built packages for the various dependencies
that I installed with something like
apt-get install xyz
put the contents in very different places than on your OSX system. The Makefile
has to change quite a bit if a package is installed in /opt vs being scattered
through /usr/local. On Windows, I don't think packages for the dependencies
even exist.

Regarding demo, I removed the non-working demo section from pde1d.m,
created an "examples" directory, and put a single file "heatCond.m" in
there.

On Tue, Jul 5, 2016 at 1:08 PM, Carlo De Falco <[hidden email]> wrote:
Bill,

On 5 Jul 2016, at 18:31, Bill Greene <[hidden email]> wrote:

> Glad to hear you were able to build on OSX. I don't currently have access to that platform.

This version of the Makefile should be relatively platform independent
as it relies on the "mkoctfile".

This is the standard approach used in all Octave packages.
Using a slightly improved version of this file you should
be able to quite easily reorganize your code into a package
using thie procedure described here:

https://www.gnu.org/software/octave/doc/interpreter/Creating-Packages.html#Creating-Packages

this is how Octave extesions are

> Regarding the demo, it looks to me like the first "demo" line is
> %!demo
> I don't see an extra space.

there shoudn't be any, the first line must be exactly

%!demo

the space should be in the other lines, as described here:

https://www.gnu.org/software/octave/doc/interpreter/Demonstration-Functions.html#Demonstration-Functions


> I haven't yet tried on linux but when I do
> demo pde1d
> on windows, I get this error:
> error: element number 1 undefined in return list
> error: called from
>     demo at line 111 column 13
>
> >Even if I fix this the demo does nothing useful, it defines a
> >bunch of functions but doesn't run anything.
>
> OK, I just took the contents of one of my test m-files and put
> a %! in front of every line. So maybe there should be
> %!heatCond
> as the last line in the demo section to call the "main" function?

typing "demo pde1d" is equivalent to copying the contents of the demo
blocks and typing "enter". If you do that, do you get any interesting result?

In the manual page for demo linked above, please read carefully the paragraph about
defining functions in a demo block:

'... Finally, because demo evaluates within a function context it is not possible to define new functions within the code. Anonymous functions make a good substitute in most instances. If function blocks must be used then the code eval (example ("function", n)) will allow Octave to see them. This has its own problems, however, as eval only evaluates one line or statement at a time. In this case the function declaration must be wrapped with "if 1 <demo stuff> endif" where "if" is on the same line as "demo" ...'

c.



Reply | Threaded
Open this post in threaded view
|

Re: 1D PDE solver for Octave

Carlo de Falco-2

On 5 Jul 2016, at 20:32, Bill Greene <[hidden email]> wrote:

> Regarding demo, I removed the non-working demo section from pde1d.m,
> created an "examples" directory, and put a single file "heatCond.m" in
> there.

I compared performance of your code to a solution of the same problem
with bim + daspk:

  >> tic, heatCond, toc
  warning: pde1d: User-defined initial conditions were changed to create a consistent solution to the equations at the initial time.
  Elapsed time is 1.26389 seconds.

  >> tic, mytest, toc
  Elapsed time is 1.05033 seconds.

on my system, though, most of the time is spent in producing the plot, so I tried
commenting out the plotting of the results in both examples, and I get:

  >> tic, heatCond, toc
  warning: pde1d: User-defined initial conditions were changed to create a consistent solution to the equations at the initial time.
  Elapsed time is 0.253297 seconds.

  >> tic, mytest, toc
  Elapsed time is 0.0677781 seconds.

This seems to confirm my previous knowledge that the mex interface is not very efficient.
If you don't want to learn Octave's C++ API, you might want to consider writing your code as
m-files. The benefits would be:

- simplify distribution and installation
- trivially solve issues with complex coefficients

I do agree that sundials is many respects superior to daspk and/or lsode so an option
to keep using that is to do so via the new implementation of ode15s instead ...

What do you think about this?

c.



mytest.m (1K) Download Attachment
Reply | Threaded
Open this post in threaded view
|

Re: 1D PDE solver for Octave

Bill Greene-3

Yes, I've been aware that this pde1d implementation is slow.
Most users probably don't care about .25 vs .07 seconds but, if the
mesh is dense, as it often needs to be with these simple linear elements,
the execution time can be many minutes; not acceptable.

I'm assuming you are probably right that the oct interface has a faster
interface than mexCallMATLAB to call user-defined m-functions. I haven't
profiled the pde1d code but I'm 99% certain that is where the time is being spent.
I'll try to find time to experiment with the oct equivalent.

Writing the code in m has some obvious attractions and if the octave
implementation of ode15s was more mature I would consider that.
You may have noticed that my code currently relies on the Sundials
ida support for banded jacobians. I don't know if octave/ode15s is
planning to support that or even how it could. What I really need for
my planned, next-version of pde1d is support for general sparse
jacobians calculated efficiently by finite difference. MATLAB/ode15s
supports this but Sundials currently does not; so, at this point, implementing
it in octave/ode15s is moot. So I'm hoping that Sundials will improve their
sparse matrix support and I expect the best way to take advantage of that
is via their c-api.

Another alternative I have been considering is adding a 'Vectorized' option.
This would work similarly to the same-named option in the MATLAB ODE functions
in that the coefficients for all mesh points would be evaluated in a single
call to the user-defined function. This would obviously require some extra
effort by the users in coding their function.

All three of these options have significant disadvantages and it is definitely
not clear to me which is the least worst.




On Tue, Jul 5, 2016 at 2:49 PM, Carlo De Falco <[hidden email]> wrote:

On 5 Jul 2016, at 20:32, Bill Greene <[hidden email]> wrote:

> Regarding demo, I removed the non-working demo section from pde1d.m,
> created an "examples" directory, and put a single file "heatCond.m" in
> there.

I compared performance of your code to a solution of the same problem
with bim + daspk:

  >> tic, heatCond, toc
  warning: pde1d: User-defined initial conditions were changed to create a consistent solution to the equations at the initial time.
  Elapsed time is 1.26389 seconds.

  >> tic, mytest, toc
  Elapsed time is 1.05033 seconds.

on my system, though, most of the time is spent in producing the plot, so I tried
commenting out the plotting of the results in both examples, and I get:

  >> tic, heatCond, toc
  warning: pde1d: User-defined initial conditions were changed to create a consistent solution to the equations at the initial time.
  Elapsed time is 0.253297 seconds.

  >> tic, mytest, toc
  Elapsed time is 0.0677781 seconds.

This seems to confirm my previous knowledge that the mex interface is not very efficient.
If you don't want to learn Octave's C++ API, you might want to consider writing your code as
m-files. The benefits would be:

- simplify distribution and installation
- trivially solve issues with complex coefficients

I do agree that sundials is many respects superior to daspk and/or lsode so an option
to keep using that is to do so via the new implementation of ode15s instead ...

What do you think about this?

c.



12