Riemann refactor for Clawpack 5.0

[Randall J LeVeque]

I'm moving this conversation to claw-dev from
https://github.com/clawpack/clawpack/issues/3
See that page first to catch up on the discussion.

Grady was wondering if putting the Riemann solver in a module and/or
using derived types causes any problems for PyClaw since f2py may have
problems with these?

When we talked in Boston there were various opinions on using derived
types. After some more thought, I'd personally be in favor of leaving
the basic calling sequence alone but adding one more parameter called
something like rpaux, perhaps, that would be auxiliary information
coming into the Riemann solver. The default definition of this type
would include the things that were passed in the common block in 4.x,
common /comxyt/ dtcom,dxcom,dycom,tcom,icom,jcom
but the user could copy this module over to an app directory and
modify it to add other things if needed.

Perhaps the user could add fields for problem parameters such as
rho,bulk,cc,zz in the acoustics example David converted, or maybe it's
better to leave those as module variables in a riemann module.

Adding other things that change from one call of the rp routine to the
next would be more complicated since the user would have to figure out
where to set these things, but for advanced applications it might be
nice to have a generic name like rpaux rather than "geometry" in order
to allow things we haven't thought of.

Grady and I were also thinking that independent of this, it would be
nice to have a geometry_module.f90 that contains utilities for dealing
with mapped grids, in particular if the user wants to recompute things
like areas and lengths in the Riemann solver rather than passing them
around in aux arrays. In 3d I think there are 13 aux variables that
would be needed just for geometry and recomputing may be a better way
to go.

We also talked about passing in xlower in addition to dx so the user
could compute x at each point in the vectorized Riemann solver, but I
think there's a problem with that in relation to the qad routine used
in amrclaw, since in this routine values are packed into the 1d ql and
qr arrays in a funny way, not just a slice of the grid, so if it's
necessary to know x or y within the Riemann solver this might not
avoid the need for putting it in the aux array. But this is a case
where maybe a derived type rpaux could be extended to include fields
xl and xr that give the x values corresponding to ql and qr that in
step2 would be computed in the obvious way using xlower and dx but in
qad would be set differently.

- Randy

[Aron Ahmadia]

> Grady was wondering if putting the Riemann solver in a module and/or
using derived types causes any problems for PyClaw since f2py may have
problems with these?

Indeed, it could be an issue for f2py, but regardless, we'll need a stable C interface, and if we have a stable C interface, we could wrap the modules with Cython instead with relatively little pain.

A

[Matthew Emmett]

If you want to bind contextual/auxiliary information at runtime, you
may consider using type(c_ptr) and c_f_pointer as shown by
https://gist.github.com/certik/1665626.

Also, if you do make a C interface, you can use ctypes (or perhaps
even CFFI) to call your Fortran routines as described by
http://fortran90.org/src/best-practices.html#interfacing-with-python
instead of have an extra Cython layer.

Matt

> --
> You received this message because you are subscribed to the Google Groups
> "claw-dev" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to claw-dev+u...@googlegroups.com.
> To post to this group, send email to claw...@googlegroups.com.
> Visit this group at http://groups.google.com/group/claw-dev?hl=en.
> For more options, visit https://groups.google.com/groups/opt_out.
>
>

[Aron Ahmadia]

Matt, that's a good point about ctypes, that's definitely the right way to approach this if we don't need any extra wrapping/types magic.

A

[Randall J LeVeque]

I took a pass at writing a Riemann module based on what I suggested,
see the comments in the attached file.

Feedback welcome!
- Randy

[Jed Brown]

On Sun, Mar 10, 2013 at 5:04 PM, Randall J LeVeque <r...@uw.edu> wrote:

I took a pass at writing a Riemann module based on what I suggested,
see the comments in the attached file.

I've talked with several people here over the years, but modules are only offering namespacing. They are still global variables, with all the implications. Most importantly, you cannot have multiple instances in the same simulation. Multiple instances are useful for multidomain problems, multigrid methods, a plethora of uncertainty quantification methods.

A while back, I helped Ondrej with this comprehensive assessment of callback mechanisms.

http://fortran90.org/src/best-practices.html#type-casting-in-callbacks

I think it would be a mistake to change interfaces without adding support for callback contexts using one of the methods above. Although Ondrej is a die-hard Fortran advocate, he also came to the conclusion that iso_c_binding offers the best way for Fortran to pass callback arguments to Fortran. Of course this is also most convenient for Python and C. The downside is that it requires Fortran 2003.

http://fortran90.org/src/best-practices.html#v-using-type-c-ptr-pointer

See also the gists at the bottom of that page.

[Randall J LeVeque]

Hi Jed,

Thanks for this link, several useful ideas on this page!

It seems like what I proposed is most similar to approach II: General
structure, in the list of ways to do callbacks, with the derived type
"rpaux" playing the role of the type named "context" there? That
seems to be one of the recommended approaches, or am I missing
something?

Is the issue the question of whether assuming the module contains
subroutines rpn2 and rpt2 with these names limits the ability to do
multidomain problems, etc? I agree that's not great. Internally
Clawpack 4.x still passes around the name of the actual Riemann solver
subroutines to use in various calling sequences, but long ago we wrote
a claw2ez wrapper that assumed they were just called rpn2 and rpt2
since that was the most frequent use, and now that's been somewhat
ossified by the fact that the Python setrun.py module doesn't specify
the name of the Riemann solver and the one actually used depend on
what files are in the Makefile.

It would be nice to allow the flexibility of having multiple Riemann
solvers used in the same run for sophisticated users who want to write
an outer program to take advantage of this, but I'd worry about making
it too complicated for the vast majority of users who don't, and not
sure we want to require Fortran 2003.

What's the best compromise?

- Randy

[Matthew Emmett]

On Sun, Mar 10, 2013 at 4:42 PM, Randall J LeVeque <r...@uw.edu> wrote:
> Hi Jed,
>
> Thanks for this link, several useful ideas on this page!
>
> It seems like what I proposed is most similar to approach II: General
> structure, in the list of ways to do callbacks, with the derived type
> "rpaux" playing the role of the type named "context" there? That
> seems to be one of the recommended approaches, or am I missing
> something?

One drawback to this approach is that the rpaux/context type will have
to be general enough to encompass all Riemann solvers, and if a user
implements a new solver they may need to modify the rpaux type.

Using the c_ptr option is extremely flexible (but perhaps at the cost
of readability).

I guess it's a matter of taste.

There was a thread last February ("workspace allocation") that
explored the same issues that are being discussed in this thread.

> Is the issue the question of whether assuming the module contains
> subroutines rpn2 and rpt2 with these names limits the ability to do
> multidomain problems, etc?

And workspace allocation too. If a multigrid algorithm tries to call
the Riemann solver multiple times with different geometries but the
workspace is static, will bad things happen?

> It would be nice to allow the flexibility of having multiple Riemann
> solvers used in the same run for sophisticated users who want to write
> an outer program to take advantage of this, but I'd worry about making
> it too complicated for the vast majority of users who don't, and not
> sure we want to require Fortran 2003.

In my experience, the GNU, Intel and Cray Fortran compilers all have
good support for the iso_c_binding module (you don't need full 2003
support, just the C binding stuff). I think this has been the case
for some years now.

Matt

[David Ketcheson]

I like Randy's proposed example.  I think we should go forward and implement this.  These discussion threads tend to fracture into a dozen different sub-issues, die off, and get resurrected a year later without any progress.

-David

[Aron Ahmadia]

I'm in favor of iso_c, which means I need to come up with an example (I'll try to put one together by the end of the month).

I'm also planning to include a proper harness for timing and exploring the various options.  Unfortunately, the traditional methods for run-time flop counting have basically been broken on Ivy Bridge (http://icl.cs.utk.edu/projects/papi/wiki/PAPITopics:SandyFlops), which is going to make it hard to count flops.

Let me know if you'd like me to run through any other ideas/examples while I'm writing code for this.

A

[Randall J LeVeque]

Aron,

I'd be interested to see an example, although I'm reluctant to make
the Riemann solvers any more complicated than necessary in terms of
what new users and students have to deal with. I'm still not sure I
understand what this would give us beyond the current proposal, and in
what contexts it is necessary, so perhaps you can illustrate the
advantages in the example.

Thanks,
- Randy

[Aron Ahmadia]

Hi Randy,

Agreed.  

Cheers,

Aron

[Kyle Mandli]

I also like Randy's example.  The issues of global module variables is fairly well contained in this case it looks like.  We could add some private declarations if this becomes a problem but I tend to agree with simplicity over flexibility in this case.

[Kyle Mandli]

After chatting with Jed a bit I see what he was trying to say about multiple instances and why using the iso_c_binding module and in particular the c_ptr type is advantageous.  The easiest use case I can think of is a problem where you may want to run multiple simulations simultaneously each with a different value of the sound speed say.  Currently one would have to write a Python script to control all of these separate runs as separate distinct processes.  This works but does limit what you can do, i.e. if you wanted to have feedback between the different realizations (stochasticity for instance).  

There's a lot more that could be said here but in the interest of brevity I will summarize as best I can in the following list:

 - I agree with Matt's statement that the c_ptr approach allows maximum flexibility at (arguably perhaps) the cost of readability.

 - The iso_c_binding support in modern Fortran compilers is robust.  My understanding is that gfortran 4.5 has a complete (as of the Fortran 2003 standard) implementation but this might be an issue, especially for Mac users who tend to grab gfortran 4.2 which has no support fot this.

 - Module globals are as bad (maybe worse) than common blocks and should be avoided.  This includes type definitions as well.  Moving the global variables in the module in Randy's example in to the definition of rpaux_type does not solve the problem (although helps a lot).

I think Aron's right that we need an example of what a c_ptr implementation of the Riemann aux data might look like although this does not alleviate the compiler versioning problem as mentioned above.

[Kyle Mandli]

Over the past few days I started playing with the iso_c_binding module a bit and have a working example:

https://github.com/mandli/classic/tree/1d-f90-conversion

There's a lot of other changes in that branch as well but the test 1d advection example should work to run and plot.  The primary point though is the Riemann solver which uses c_ptr to transfer arbitrary data to through to the Riemann solver and a geometry type which is passed which contains x, t, dx, dt.  In this particular case a point-wise Riemann solver is being used but this can be swapped out easily.  Anyone have any thoughts?  I have a really simple example of the usage of what is needed to do this if anyone is interested in a simpler example.

Kyle

[Aron Ahmadia]

Hi Kyle,

Thanks for this, and sorry for the delay in responding.  At first I was looking at your last commit, but it appears that you've got quite a bit of other refactoring against the classical interfaces in the other commits, so it might be easier if I spend some time getting a clearer diff together.  I'll try to summarize some of your changes/rationale here.

https://github.com/mandli/classic/commit/fbcd48ec1b724b6d561cb9eba36c57540473210d

One a file-by-file basis let me try to summarize:

New file: src/1d/geometry_module.f90

In this file, you define the geometry of the grids that this Riemann solver will operate on:

+    type geometry_type

+        real(kind=dp) :: x, dx, t, dt

+    end type geometry_type

I believe x here is some sort of physical indication for the lower-left-hand corner of the grids, t is the base time, and dx/dt are obvious :)

It makes sense that if we modify the Riemann routines, we're going to have to modify the interfaces :)

Modified file: src/1d/hyperbolic_step.f90

        do i=2-solver%num_ghost,solution%num_cells(1)+solver%num_ghost

            geometry%x = solution%lower(1) + (i-0.5_dp) * solution%dx(1)

            geometry%dx = solution%dx(1)

            geometry%t = solution%t

            geometry%dt = solver%dt

            call solver%rp1(solution%num_eqn,               &

                            solution%num_aux,               &

                            solver%num_waves,               &

                            solver%rp_data,                 &

                            geometry,                       &

                            q(:,i-1), q(:,i),               &

                            aux(:,i-1), aux(:,i),           &

                            wave(:,:,i),                    &

                            s(:,i),                         &

                            amdq(:,i),                      &

                            apdq(:,i))

        end do

Here we see the geometry module being populated, then pushed into a solver call.   You could have also done one step better and lifted geometry out of the loop:

        ! lift geometry assignment out of inner loop

            geometry%x = solution%lower(1) + (i-0.5_dp) * solution%dx(1)

            geometry%dx = solution%dx(1)

            geometry%t = solution%t

            geometry%dt = solver%dt

        do i=2-solver%num_ghost,solution%num_cells(1)+solver%num_ghost

            call solver%rp1(solution%num_eqn,               &

                   ...

        end do

It would, of course, have been nice if we could have pushed that data directly into an object or some other static location specific to the function, but in Fortran, we're limited by the language features we're willing to require from our compilers.  Since we're not going to require type-bound procedures (Fortran 2003), we're really only left with a ISO C compatible (also 2003, but apparently better-supported) derived type pointer that acts as a bucket for any data that does not go into the mandatory interface.

I've got more thoughts (and more code reviewing to do), but I'd like to spend some time getting some test harness code together in Riemann so everybody can play with the different approaches.  

More soon,

Aron

[Kyle Mandli]

Thanks Aron for taking a look, this set of changes are definitely a departure from what was previously implemented and the branch name is not really indicative (refactor of the code) of what is going on.  The major changes here are meant to increase readability of the code and encapsulation of data (no module variables if possible). 

One additional aspect of this that we should consider is f2py's support for the iso_c_binding module.  I have not spent a lot of time looking into this but was not able to use f2py to compile code that used the c_ptr type.  f2py defaults to using void * when it cannot figure out what is going on so it seems like it might be possible but it might require an interface file.

Since this was implemented using a point-wise Riemann solver this is the natural way to do this but you are probably right that this could be done in a vectorized form outside of the loop and only indexed into.

[Aron Ahmadia]

All,

I'm going to recap what's happened in the last year on this.

I'll use static to refer to the existing style of interfaces to Riemann solvers, and object-based to refer to the interfaces that Randy, Jed, Matt, and Kyle have been proposing.  At this point, we have decided on a basic refactor to the static Riemann solver interfaces, and I'm trying to come up with an object-based code solution that will keep everybody happy.  Just to reiterate the existing discussion/options one more time:  

* The existing static solvers provide a single implementation of the Riemann solver code as a subroutine using only Fortran native types, and access all solver-specific data through function parameters such as the aux arrays, which are of arbitrary dimension, and common blocks, which are static data structures shared across the program's address space.  We agreed to a static solver refactor at SIAM CSE, but we didn't do a clear job of identifying how we would implement an object-based approach.  I think several people tried to raise this, and I apologize for the part I played in skipping this aspect of the discussion.

David did a great job of putting together an example of the basic static refactor we agreed on, using module data instead of a common block: https://github.com/ketch/riemann/blob/riemann_modules/src/rp1_acoustics.f90

* object-based solvers provide a single implementation of the Riemann solver code as a subroutine using Fortran native types and an extensible generic type, with no common blocks.  Solver-specific data is passed in through the generic type, enabling multiphysics solvers and other situations where having multiple solvers is needed.  As I mentioned earlier, Jed and Matt both have use cases for this in their existing applications.  We've bounced around the idea of refactoring to an object-based approach, and Randy, Kyle, and Matt (with constructive feedback from Jed) have provided examples of how this could be done:

    + Randy's example: https://gist.github.com/ahmadia/5382059

    + Kyle's example: https://github.com/mandli/classic/blob/fbcd48ec1b724b6d561cb9eba36c57540473210d/tests/advection/rp1_advection.f90

    + Matt's example: https://github.com/memmett/pyclaw/blob/memmett-ctypes-pdata/apps/acoustics_1d_homogeneous/rp1_acoustics.f90

As I've been re-learning, Fortran has fairly weak support for the concept of objects, with derived types being the closest thing to proper support for what we need for native Fortran generic data structures that represent Riemann solver data.  In all three code examples, a derived type is created for handling the solver-specific data, and in Kyle's/Matt's examples, you can see how this code is further wrapped into a C type.  Matt's example (from February 2012!) is particularly nice, as the complete tree at: https://github.com/memmett/pyclaw/tree/memmett-ctypes-pdata/apps/acoustics_1d_homogeneous gives us a reasonable path forward for integration into PyClaw.  The unfortunate aspect of this is that if we switch back to ctypes, we have to put Makefiles *back* into the code base.  One option is to integrate them into the build process in the same way Lisandro does in petsc4py and mpi4py.  Another option would be to have a completely separate Clawpack 5 library  

There were some concerns raised about performance and what changes would need to be made above the Riemann solvers to support the new interfaces.  I've created some test code in PyClaw to hopefully address both of these:

https://github.com/ahmadia/clawpack/tree/riemann_shootout

I'm planning on integrating Matt Emmett's ISO C approach, since it's the most complete, but if there are salient differences between it and Kyle or Randy's code, please let me know so I can create another variant.  I'll post again when I've got Matt's code integrated, but if you're interested in playing with a simple comparison between native Fortran and Python versions of the 1d advection code, you can install the riemann_shootout code with the following commands:

git clone g...@github.com:ahmadia/clawpack.git

cd clawpack

git checkout -t origin/riemann_shootout

pip install -e .

cd pyclaw

# verify you've got a working pyclaw build

nosetests 

cd apps/riemann_shootout

python riemann_compare.py

Cheers,

Aron

[Randall J LeVeque]

Hi Aron,

Thanks for keeping the discussion alive. When Kyle was here a couple
weeks ago we talked more about this and I'm coming around to the view
that iso_c_bindings might be the best way to go. I also think it's
probably a good idea to have a separate geometry parameter that can
always be assumed to contain what might be needed in a Riemann solver
for cell geometry, e.g. cell center and dx,dy in computational space
at least, from which the corners/edges could be computed if needed for
a grid mapping. These are easy enough to set in the calling routine,
e.g. flux2, and used often enough to be useful to make the default.
In this regard it should be remembered that in amrclaw the Riemann
solver is called also from subroutine qad (for conservation fix up)
with ql and qr not being a set of adjacent states along a single row
of the grid. Instead these are filled with values along the edge of a
fine grid patch on one side and from adjacent cells in the coarser
parent patch on the other side.

So I think it's a good idea to have a geometry argument that would be
of a fixed derived type and also something like Kyle's rp_data or
Matt's context that would be an iso_c pointer to whatever else the
user might want to pass through to the solver. I believe this would
make it relatively easy for most users - if they are not doing
something advanced where they need rp_data then they wouldn't have to
learn about iso_c_bindings. Correct me if I'm wrong, but I would
guess if rp_data isn't used in the Riemann solver then it could just
be a dummy argument of some arbitrary type and the user wouldn't have
to "use" the iso_c_binding module or worry about having an appropriate
compiler? Or does it have to be declared properly in library routines
like flux2 that it is passed through, even if unused?

So I think I like Kyle's proposal, but I haven't had a chance to
experiment with it for myself yet.

I'd be happy to hear feedback from others on the various options now
on the table.


- Randy

[Aron Ahmadia]

All,

As of the following commits:

ahmadia/clawpack:

48711385e27db23e85c2c0e16047c31334521afe

clawpack/pyclaw:

5de2750128e9c8b1ada2f77640c39584437b77b4

There is a new demonstration of the ctypes functionality calling into unwrapped Fortran library code using ISO C bindings.  Now you can test an advection kernel being run through native Python, the standard "Fortran" interface to PyClaw, and a proposed ISO C variant.  You should also have a look in:

clawpack/pyclaw/apps/riemann_shootout/next

and

clawpack/pyclaw/apps/riemann_shootout/next/iso_c

for a feel on the changes to the code.  I resisted the urge to make many modifications outside of the bindings.  I haven't figured out yet which variables belong in the derived type vs. the calling interface, if somebody would like to see another particular interface tried out, I'm happy to try and code it up.

I mostly followed Matt Emmett's approach, and I was able to implement the new interfaces as a new 1D classic solver subclass.  I don't think it will be very hard to extend this approach beyond the prototype implementation.

I will post a more detailed email later describing how things work, but it would be helpful if others could download and test on their own development machines:

# in a new directory

git clone g...@github.com:ahmadia/clawpack.git

cd clawpack

git checkout -t origin/riemann_shootout

# this will be less dangerous if you're working from a virtualenv, but do what works for you

pip install -e . # or some variant of python setup.py build/install

cd pyclaw

# verify you've got a working pyclaw build

nosetests 

cd apps/riemann_shootout

# build the ISO C bindings

cd next

# you may need to adjust your Makefile, use f2py -c --help-fcompiler for guidance on adjusting your compile flags

make

cd ..

python riemann_compare.py

One current issue is that my modifications have introduced a 30-50% slowdown over the wrapped f2py approach.  I haven't had a chance to look deeply into the cause, but it's on the list of to-dos.

A

[Jed Brown]

Aron Ahmadia <ar...@ahmadia.net> writes:

> There is a new demonstration of the ctypes functionality calling into
> unwrapped Fortran library code using ISO C bindings. Now you can test an
> advection kernel being run through native Python, the standard "Fortran"
> interface to PyClaw, and a proposed ISO C variant.

Aron, why does this use iso_c types everywhere?

If I was writing this primarily for Fortran users, I would use Fortran
types everywhere except when _library_ code needs to pass the context.
When user code gets the context (as a c_ptr), I would immediately
promote it to its Fortran type, and I would only reduce it to
c_ptr(param) when passing it as an argument to a library function.

Now if your objective is to export an ABI that is usable without shims
from C, Python, and Fortran, what you have is a reasonable approach, but
I think it's unnatural as a Fortran API and I think that trying to mix C
and Fortran without shims tends to produce crummy interfaces on both
sides. As someone that prefers programming in C, I won't object to
making it more convenient for me to implement Riemann solvers in C, so
if you're happy with the somewhat "foreign" Fortran interface, I won't
try to stop you.

[Aron Ahmadia]

Hi Jed,

Thanks for taking a look at the code and sorry for the delay in response, I was off facing 50 mph winds and hail on top of Cadair Idris this weekend :)

> Aron, why does this use iso_c types everywhere? ... Now if your objective is to export an ABI that is usable without shims from C, Python, and Fortran, what you have is a reasonable approach, but I think it's unnatural as a Fortran API and I think that trying to mix C and Fortran without shims tends to produce crummy interfaces on both sides

Yes, interoperability was my objective.  Because I need interoperable procedures at the outer kernel level and the Riemann kernel level, I modified the interfaces to both the outer solver (step1.f90) and the Riemann kernel (rp1_advection.f90)

> If I was writing this primarily for Fortran users, I would use Fortran

types everywhere except when _library_ code needs to pass the context.
When user code gets the context (as a c_ptr), I would immediately
promote it to its Fortran type, and I would only reduce it to
c_ptr(param) when passing it as an argument to a library function.

The "library" code in your statement is the outer solver routine that calls the Riemann solver kernel (step1.f90 and friends), and your 'user' code is the Riemann kernel itself (rp1_advection.f90 and other Riemann kernels).  Just to make this more precise, I'm going to use "outer solver", "Riemann kernel", and a third type of code, "driver", which refers to code that sets up the problem domain and solver routines before calling them.  

My understanding is that unlike PETSc, there is pretty equal development in all three types of code, and a new user is equally likely to be working with any of the three code bases, depending on whether they are focusing on existing applications, time integrators, or solving new applications.  Of the three sections of code, only the "outer solver" deals with the context as an opaque type.  If you look at step1.f90, you'll notice that it leaves context as an opaque type throughout the code.  One aspect of the code you may be taking issue with is my choice of passing a C function pointer in to represent the Riemann kernel instead of the native Fortran function pointer? I agree that I could have written the core routine to be native Fortran, but since we are already incorporating ISO C bindings in this routine's interface, we might as well bind the routine to make it interoperable instead of writing a redundant shim function.

In the Riemann kernel rp1_advection.f90, I immediately move from the opaque c_ptr to a Fortran-native derived type and use it there.  The same interface and module is used in the driver routine to set up and destroy the context object (the new and delete methods).  I think this is what you're suggesting, so I'm confused as to where exactly your criticism is being directed at.  I suspect I'm still missing something here, so I'd appreciate clarification when you've got a few minutes.

Cheers,

Aron

[Aron Ahmadia]

Ondrej's response here on StackExchange is relevant to this discussion, in particular because he has a lot of other advice on modernizing Fortran code that goes beyond my own expertise: 

http://scicomp.stackexchange.com/questions/6961/f2py-with-allocatable-and-assumed-shape-arrays

A

[Aron Ahmadia]

Thanks for your input Ondrej!

I've been cropping this thread in my email responses, but you can see it in its full glory here: https://groups.google.com/forum/?fromgroups=#!topic/claw-dev/7G67NxwzxsE

I tried to recap the last year or so of activity on this here: https://groups.google.com/d/msg/claw-dev/7G67NxwzxsE/rOJ4UtgMvqIJ

The interface changes I'm proposing are here: https://groups.google.com/d/msg/claw-dev/7G67NxwzxsE/rOJ4UtgMvqIJ

The Clawpack library is not very big, and I was hoping I could get away with just standardizing the interfaces on ISO C, which would simplify my job, but if both you and Jed are strongly suggesting stub interfaces, it looks like I'll need to write them over a native Fortran API using derived types.  

I think there is a large contingent of Fortran developers on claw-dev who will be glad to hear that you advocate working with Fortran instead of a mixed-language approach.  My understanding is that there is still some hesitation to embrace a lot of the modern Fortran programming techniques, mostly out of concern for portability and performance.  I'm glad resources such as your guide exist, because you spend a fair amount of time addressing these concerns.

That said, we really like Python :)

Cheers,

Aron

On Thu, Apr 25, 2013 at 12:18 AM, Ondrej Certik <ondrej...@gmail.com> wrote:

Hi,

since Aron CCed me, I decided to write my personal experience with Fortran.
Maybe you might find it useful. But I am not suggesting anything for the clawpack project, I am only posting it here for inspiration:

On Monday, April 22, 2013 4:25:03 AM UTC-6, Aron Ahmadia wrote:
> Hi Jed,
>
>
> Thanks for taking a look at the code and sorry for the delay in response, I was off facing 50 mph winds and hail on top of Cadair Idris this weekend :)
>
>
>
> > Aron, why does this use iso_c types everywhere? ... Now if your objective is to export an ABI that is usable without shims from C, Python, and Fortran, what you have is a reasonable approach, but I think it's unnatural as a Fortran API and I think that trying to mix C and Fortran without shims tends to produce crummy interfaces on both sides
>
>
>
> Yes, interoperability was my objective.  Because I need interoperable procedures at the outer kernel level and the Riemann kernel level, I modified the interfaces to both the outer solver (step1.f90) and the Riemann kernel (rp1_advection.f90)

In general, I think the best approach is to code in modern Fortran (that is, no C types), and use everything that modern Fortran offers (as long as it is reasonably supported by compilers), in particular array support, allocatable arrays etc., all of which you loose if you stick to iso_c_binding, which only has explicit shape arrays.


But it depends how complicated your program is. If it is just one subroutine, then your approach might make sense. But if it is more subroutines, calling each other (that is typical in my programs/libraries), then I would recommend to write a simple C wrapper for the API that you want to export.

>
>
>
> > If I was writing this primarily for Fortran users, I would use Fortran
> types everywhere except when _library_ code needs to pass the context.
>
> When user code gets the context (as a c_ptr), I would immediately
> promote it to its Fortran type, and I would only reduce it to
>
> c_ptr(param) when passing it as an argument to a library function.

Exactly.

>
>
>
> The "library" code in your statement is the outer solver routine that calls the Riemann solver kernel (step1.f90 and friends), and your 'user' code is the Riemann kernel itself (rp1_advection.f90 and other Riemann kernels).  Just to make this more precise, I'm going to use "outer solver", "Riemann kernel", and a third type of code, "driver", which refers to code that sets up the problem domain and solver routines before calling them.  
>
>
>
> My understanding is that unlike PETSc, there is pretty equal development in all three types of code, and a new user is equally likely to be working with any of the three code bases, depending on whether they are focusing on existing applications, time integrators, or solving new applications.  Of the three sections of code, only the "outer solver" deals with the context as an opaque type.  If you look at step1.f90, you'll notice that it leaves context as an opaque type throughout the code.  One aspect of the code you may be taking issue with is my choice of passing a C function pointer in to represent the Riemann kernel instead of the native Fortran function pointer? I agree that I could have written the core routine to be native Fortran, but since we are already incorporating ISO C bindings in this routine's interface, we might as well bind the routine to make it interoperable instead of writing a redundant shim function.
>
>
>
> In the Riemann kernel rp1_advection.f90, I immediately move from the opaque c_ptr to a Fortran-native derived type and use it there.  The same interface and module is used in the driver routine to set up and destroy the context object (the new and delete methods).  I think this is what you're suggesting, so I'm confused as to where exactly your criticism is being directed at.  I suspect I'm still missing something here, so I'd appreciate clarification when you've got a few minutes.

Is this the code that we are talking about?

https://github.com/clawpack/riemann/tree/master/src
https://github.com/clawpack/pyclaw/blob/master/src/pyclaw/classic/step1.f90

-------------------------


I personally have been using Python for probably 13 years now, so I am reasonably proficient in it.
I used to have the main thing in Python and call little C kernels from it.
Then I learned modern Fortran, about 2 years ago. Since then I actually just use modern Fortran only, because to me it feels like Python/NumPy, just with types (types and compiler can really help!), and much much faster, for free.
And there is high value in using just one language, rather than a mix.

I wrote this Python/Fortran side by side comparison page:

http://fortran90.org/src/rosetta.html


Here is my Fortran code for inspiration:

https://github.com/certik/hfsolver/tree/master/src

and here is another code (F95 for maximum compatibility) that we published in Computer Physics Communications recently (article is cited in the README):

https://github.com/certik/dftatom

In dftatom I do have optional Python wrappers, but only for the high level functionality.


My experience from the last 2 years is this (this is just my own opinion and I know that some people will not agree with me):

* it's really simple to write very fast code (usually faster than equivalent C++ code, see e.g. my recent test at [1])

* modern Fortran code looks similar to the mathematics it solves (arrays, complex numbers, ...), that really helps

* Debugging and profiling Python+Fortran mix is not worth my time, it's much faster to just use Fortran for everything. Both the result and my development time is faster that way. In particular, the dftatom used to be 100x slower than it is now. We needed it to be as fast as an old legacy optimized f77 codes (see the article for details). It was really tough --- I even suspect that f77 might sometimes be a little faster than f90, but only slightly. In particular, I had to give up Python and just put everything in Fortran to get the job done in my time frame. Things became simpler, faster and more maintainable.

* For any nontrivial numerical code that can be well expressed using arrays (e.g. almost all algorithms that I use), if I need to squeeze the maximum performance, then I would definitely choose Fortran over C or C++.
One such example is this competition that we did some time ago [2].
The goal was to solve some mathematical problem the fastest.
First we all use different algorithms, so it wasn't clear which language is the fastest. But then Jed figured out the fastest algorithm. He coded it in Octave and Julia. Aron coded it in Cython/NumPy. I coded it directly in modern Fortran.
All three versions are very readable and pretty much equivalent, as all Octave, Julia, NumPy and (modern) Fortran have pretty much the same expressiveness.
Nobody tried C or C++, as probably it would take much longer to write.

And the Fortran version beats all other approaches by an order of magnitude.
And that's pretty much my whole point --- if I need to get the fastest code, and I have little time (just like everybody else), then Fortran is the way to go. Also it shows why Python+Cython or NumPy is many times not competitive in my experience if I need raw speed for nontrivial programs.

* I like to use Python for experiments and plotting. I use ipython notebook for that.

* Easy to compile, runs on clusters etc. However, as a Python enthusiast, I do want to help with the Python compilation problem, I do believe it doesn't have to be that hard as it is now to compile mixed Python/C/Fortran code.


Jed and I had a long debate about modern Fortran about a year ago I think. Let's say that he didn't share my enthusiasm for modern Fortran. :) But the above is my honest experience.

Again, I just wanted to point my own experience, and I am in no way suggesting anything.

Ondrej

[1] http://scicomp.stackexchange.com/questions/5477/evaluate-the-sum/6831#6831

[2] https://plus.google.com/u/0/116518787475147930287/posts/jiULACjiGnW

[Jed Brown]

Aron Ahmadia <ar...@ahmadia.net> writes:

> The "library" code in your statement is the outer solver routine that calls
> the Riemann solver kernel (step1.f90 and friends), and your 'user' code is
> the Riemann kernel itself (rp1_advection.f90 and other Riemann kernels).
> Just to make this more precise, I'm going to use "outer solver", "Riemann
> kernel", and a third type of code, "driver", which refers to code that sets
> up the problem domain and solver routines before calling them.

I'm happy to include the Riemann kernels as "library code", they are
just a different library component from the outer solver. In
particular, Riemann solvers implement an interface specified by the
outer solver, but neither has any further knowledge of the other.

> My understanding is that unlike PETSc, there is pretty equal development in
> all three types of code, and a new user is equally likely to be working
> with any of the three code bases, depending on whether they are focusing on
> existing applications, time integrators, or solving new applications. Of
> the three sections of code, only the "outer solver" deals with the context
> as an opaque type. If you look at step1.f90, you'll notice that it leaves
> context as an opaque type throughout the code. One aspect of the code you
> may be taking issue with is my choice of passing a C function pointer in to
> represent the Riemann kernel instead of the native Fortran function
> pointer? I agree that I could have written the core routine to be native
> Fortran, but since we are already incorporating ISO C bindings in this
> routine's interface, we might as well bind the routine to make it
> interoperable instead of writing a redundant shim function.

That gives up type safety and makes the interface more awkward. The
"outer solver" needs the type to be opaque because it is used
generically. But the driver (as you write them, since you don't have a
run-time configuration system -- PETSc's does what Martin Fowler calls
"service locator") knows which Riemann solver they are interacting with,
and in fact, they are responsible for configuring that Riemann solver,
so it's better for them to use the typed interface rather than passing
around a c_ptr. This is the reason PETSc uses strong typedefs for
objects that you interact with and void* only for generic arguments.

Now if you want people to be able to write generic drivers (e.g., one
driver that can be used for shallow water, Euler, MHD, etc; maybe
configured using a GUI), then the Riemann solver interface would need to
become richer. In particular, you would need a generic Riemann solver
interface that can be used by the driver. This would likely be a
thicker interface than that used by the outer solver and may justify
using inheritance/virtual functions.

In PETSc, we prefer passing (function_ptr, context) when we don't want
to dictate the method of construction and lifecycle, but create a real
object (with delegator-style dispatch) when we want a runtime-configured
object that can be created generically and is reference counted. Since
a generic driver needs uniform lifecycle management, passing
(function_ptr, context) would be more cumbersome.

[Aron Ahmadia]

Hi Jed,

Thanks for your points.

> I'm happy to include the Riemann kernels as "library code", they are

just a different library component from the outer solver.  In
particular, Riemann solvers implement an interface specified by the
outer solver, but neither has any further knowledge of the other.

Yes, spot on.

> That gives up type safety and makes the interface more awkward.  The
"outer solver" needs the type to be opaque because it is used
generically.  But the driver (as you write them, since you don't have a
run-time configuration system -- PETSc's does what Martin Fowler calls
"service locator") knows which Riemann solver they are interacting with,
and in fact, they are responsible for configuring that Riemann solver,
so it's better for them to use the typed interface rather than passing
around a c_ptr.  This is the reason PETSc uses strong typedefs for
objects that you interact with and void* only for generic arguments.

There are two forces acting in opposition here.

(A) A generic/extensible type that can freely be passed between the driver (full type), outer solver (abstract), and the inner Riemann kernel (full type).

(B) ISO C bindings for the functions and their data.  

What I'm still trying to wrap my head around is how to bind (B) to (A), i.e., expose a C API for all three sections of the code without deep copies.  We can't directly provide C bindings to an extended type using ISO C Bindings, so my proposed solution was to provide specific constructors/destructors for the Riemann kernel-specific context from the Riemann kernel (see the new/delete methods in rp1_advection.f90).  If you or Ondrej know of a better (type-safe?) way to provide an extensible Fortran type that can be passed either opaquely or fully typed interoperably between C and Fortran, that would be great!

> Now if you want people to be able to write generic drivers (e.g., one

driver that can be used for shallow water, Euler, MHD, etc; maybe
configured using a GUI),

Let's stay away from that discussion for now.

> In PETSc, we prefer passing (function_ptr, context) when we don't want

to dictate the method of construction and lifecycle, but create a real
object (with delegator-style dispatch) when we want a runtime-configured
object that can be created generically and is reference counted.  Since
a generic driver needs uniform lifecycle management, passing
(function_ptr, context) would be more cumbersome.

I understand in C that we could use a simple struct for the stub type, which provides a strong type for the user to interact with, and an opaque pointer within that holds all of the extended data.  

Are you suggesting that we create something similar for the Riemann kernel context?  

! abstract type

type, bind (C) :: rp1_data

    type (c_ptr) :: context

end type 

! implemented parameters

type :: advection_context 

    real(dp) :: u

end type

subroutine advection(...., rp_data) 

...

type(rp1_data), intent(in), value :: rp_data

...

type(advection_context), pointer :: advection_param

call c_f_pointer(rp_data%context, advection_param)

! now use advection_param%u wherever you like in kernel

I think this would make the Fortran driver and Riemann kernel interfaces slightly safer, at the expense of some extra code.  My big gripe is that it will require building/generating/maintaining C interfaces separately from the Fortran interfaces, and I would still rather maintain less code by sticking to a single interface which can be directly implemented in either C or Fortran.  If everybody else thinks it's worth doing, though, I'm not going to fight it, and I'm happy to create another version that uses this strategy.

A

[Aron Ahmadia]

Also, could somebody who knows more about Fortran than me respond to Randy's earlier question more effectively?

"Correct me if I'm wrong, but I would

guess if rp_data isn't used in the Riemann solver then it could just
be a dummy argument of some arbitrary type and the user wouldn't have
to "use" the iso_c_binding module or worry about having an appropriate
compiler?  Or does it have to be declared properly in library routines
like flux2 that it is passed through, even if unused?"

I think this gets to the heart of the question on binding an extensible type in Fortran, and as far as I know, the answer is: You will need a freely available compiler released in the last 6 years to be able to stably work with code that uses ISO C bindings.  Once the Riemann code base adds the ISO C bindings to wrap the kernel context, the code will not compile on older compilers, and there's no way to separate the C bindings out from the wrapped code.  

 I don't think the context needs to be appropriately declared in outer solver like flux2/step2, but my suggestion would be to convert all of the inner Riemann kernels (and eventually the outer solver routines) to export a module interface (in addition to the C interface) so that if somebody wants to write type-checked code, it is available to them.

A

[Kyle Mandli]

I think you could get away with not typing the rp_data correctly if it
was not explicitly defined which is how work arrays were repurposed
into different shaped arrays. Unfortunately if we put the Riemann
solvers into a module this automatically makes them explicitly defined
and anytime someone calls them from somewhere else the compiler will
check the arguments which means you would have to include the
iso_c_binding module.

[Jed Brown]

Aron Ahmadia <ar...@ahmadia.net> writes:

> You will need a freely available compiler released in the last 6 years
> to be able to stably work with code that uses ISO C bindings. Once
> the Riemann code base adds the ISO C bindings to wrap the kernel
> context, the code will not compile on older compilers, and there's no
> way to separate the C bindings out from the wrapped code.

CPP macros. ;-)

Dummy parameters are never dereferenced so you only need to make the
type checker accept it, readability be damned.

[Kyle Mandli]

Having had a chance to look at what Aron has, I feel somewhat inclined to say that we may be giving up too much in terms of readability to maintain full C compatibility directly without separate shims.  In terms of the slow down, is this due to copies?  I have very little knowledge on how all this works in terms of the C-type compatibilities with arrays.

As for CPP macros, I shudder at the thought, makes me feel like the old man on the porch yelling "get ou' of my code yah hooligans!  Kids these days with all of their new-fangled macros!" :P

[Jed Brown]

Kyle Mandli <kyle....@gmail.com> writes:

> As for CPP macros, I shudder at the thought, makes me feel like the old man
> on the porch yelling "get ou' of my code yah hooligans! Kids these days
> with all of their new-fangled macros!" :P

I was not being serious, if there was any question...

[Jed Brown]

Aron Ahmadia <ar...@ahmadia.net> writes:

> I understand in C that we could use a simple struct for the stub type,
> which provides a strong type for the user to interact with, and an opaque
> pointer within that holds all of the extended data.

Much more natural in C is:

typedef struct Riemann_private *Riemann;

Then we interact with 'Riemann' everywhere in the code, and only the
constructor allocates the struct Riemann_private. Since the user never
sees the contents of Riemann_private, we can do anything inside there
without changing the ABI.

> Are you suggesting that we create something similar for the Riemann kernel
> context?
>
> ! abstract type
> type, bind (C) :: rp1_data
> type (c_ptr) :: context
> end type
>
> ! implemented parameters
> type :: advection_context
> real(dp) :: u
> end type
>
> subroutine advection(...., rp_data)
> ...
>
> type(rp1_data), intent(in), value :: rp_data
> ...
>
> type(advection_context), pointer :: advection_param
> call c_f_pointer(rp_data%context, advection_param)

I don't see any value in this indirection. Since the driver knows
exactly which Riemann kernel it's using, we can make it simple:

type(advection_context), target :: adv
adv%u = 2.18

If we were to call directly, it would look like:

call advection_rp1(..., c_loc(adv))

but normally the discretization ("outer solver") is calling this. In
that case, the driver would call something like:

solver_set_riemann(solver, advection_rp1, c_loc(adv))

You can use more objecty syntax here if you prefer. The kernel itself
looks like:

subroutine advection_rp1(..., c_ctx)
type(c_ptr), intent(in) :: c_ctx
type(advection_context), pointer :: adv
call c_f_pointer(c_ctx, adv)
...
flux(i) = adv%u * q(i)

[Aron Ahmadia]

Okay, I understand.

The approach you are advocating is an opaque C pointer resolving to a (potentially unbound) Fortran type, as opposed to an opaque C pointer resolving to a ISO C bound Fortran type.

This indeed hides almost every aspect of the ISO C bindings from the core Fortran routines.  Additionally, since the context itself is simply a C pointer, if we were to implement C-only (or another language) routines in the future, we could provide a context that makes sense for that language.

Thanks, this is indeed far more flexible than the approach I was proposing.  I'll modify the example code and start a new thread for people to look at.

A 

[Randall J LeVeque]

Hi Ondrej -- Sorry it's taken so long to respond, this quarter has
been crazy. But thanks for the input and the plug for Fortran. Many
of us will certainly continue to use it, and the pages you've put
together on best practices and the Rosetta are great resources as
well. I forwarded these links to the class I'm teaching, where some
students are faced with learning both Python and Fortran
simultaneously.

All -- Regarding Riemann solvers in Clawpack, I've talked offline with
several others on claw-dev recently about how to move forward. I
suggest that in order to have nonzero probability of releasing 5.0
this summer we should stick with the Riemann solvers currently in the
riemann repository for now. So I suggest that if others want to
develop new versions for PyClaw, e.g. with iso_c_bindings and/or use
of derived types and module variables, that should be done in a new
subdirectory of the riemann repository (and perhaps move the current
version into a named subdirectory so the two or more versions are
parallel).

There are enough other things to get working for 5.0 without trying to
change all the Riemann solvers at the same time, and I'd like to get
5.0 out so that various other new features and bug fixes can be
applied to that rather than to 4.x. Then once that is stable and
producing results consistent with 4.x we can think about a potentially
massive change to the Riemann solvers for the next major release.

I've chatted with Aron, David, Kyle, Marsha, and a few others about
this and it seems like the best way to go, but if others want to chime
in please do.

I've also been slow to get the claw-dev hackathon organized, but we're
still planning this for the week of July 22 in Seattle. The tentative
plan is for some informal talks the first few mornings on topics of
mutual interest to help us move forward, but most of the time devoted
to hacking. Details to follow.

- Randy