2
$\begingroup$

I have been writing a software package in Fortran for solution of the Vlasov-Poisson system in 2D2V. I want this software to be useful beyond its current application (e.g. systems with different number of position or velocity dimensions, change in numerical methods used, solution of full Maxwell system, etc.), which led me to using abstract data types (ADTs) for some of the underlying structures to improve generality and reusability.

In my current discretization, I am using stretched Cartesian grids in space and Hermite basis functions in velocity. Here's a sketch of my current implementation:

type cartesian_grid
    ! dimension of grid
    int :: N 
  
    ! physical coordinates
    double precision, pointer :: x(:) 

    ! metric
    double precision, pointer :: xi_x(:) 

    ! step size of computational coordinates
    double precision :: delta_xi  
end type
type hermite_grid
    ! length of coordinate array
    int :: N  

    ! coordinate array
    double precision, pointer :: v(:)
  
    ! Moments of basis functions - used to compute moments of f.
    double precision, pointer :: I0(:), I1(:)  

    ! Scaling factor of Hermite basis functions
    double precision :: ux 
end type

For my time integration, I need to compute the RHS of the Vlasov equation: $$ \partial_t f=-v\cdot\nabla_rf-\frac{q}{m}E\cdot\nabla_vf $$ Currently, this is done with a subroutine called like so:

call vlapo_dfdt( f, cartesian_grid1, cartesian_grid2, hermite_grid1, hermit_grid2, dfdt)

The last component of this is an explicit time integrator, which integrates $f$ from $t$ to $t+h$. I have written number of time integrators (RK4, RK Feldberg, AB2, etc.). These are currently called using the following interface:

call time_integrator( t, h, f, cartesian_grid1, cartesian_grid2, hermite_grid1, hermite_grid2, vlapo_dfdt)

Here is where I run into trouble. I'd like to generalize this so that the grid representations are not fixed. This way I can change the velocity grids to (say) a Cartesian representation without having to write another time integration routine. The trouble is that my current integrator expects a type(hermite_grid) and won't accept a type(cartesian_grid) for the velocity grids.

I want to generalize this so that I integrate $f$ from $t$ to $t+h$ as follows:

call time_integrator( t, h, f, space_grid1, space_grid2, vel_grid1, vel_grid2, deriv_function)

where space_grid1, space_grid2, vel_grid1, vel_grid2 can be either Hermite or Cartesian grids and deriv_function is a procedural variable with the same interface as vlapo_dfdt(). This way I could change the implementation of the derivative without needing to re-write the time-steppers.

It seems like the problem is that my ADTs are not abstract enough but how can I further generalize them while still passing the data needed to compute the derivative into my deriv_function()? For cartesian coords I need the metric and computational coord step size. For Hermite coords I need the moments of my basis functions and scaling factor. For other representations I'd need other data altogether.

What is the best way to accomplish this in modern Fortran? This seems like such a ubiquitous problem in the computational sciences that I assume there's a standard (or best practice) way to do this.

$\endgroup$
2
$\begingroup$

It sounds like what you want is an abstract type -- a common grid type that both Cartesian and Hermite grids inherit from. This abstract grid type doesn't actually do anything by itself, but rather it defines an interface that callers can count on the actual implementations, namely Cartesian and Hermite grids, to provide. Ideally, you'd be able to code higher-level routines, like time integration, in such a way that the details of the lower-level entities like the grid are not important. This is an example of what's called the template method pattern.

Fortran 2003 and later are fairly standard single-dispatch on type object-oriented language, similar to Java before generics. The mechanics of how you actually implement and use abstract base classes are described here with some code samples. If you want to read more, this book by Damian Rouson is quite good. I also recommend reading through fortran90.org which has a lot of great information.

Leaving aside the specifics of exactly how you implement this in Fortran, the Design Patterns book is a must-read if you do any object-oriented programming regardless of what language you're working in. The template method is one of the patterns they describe but there are many others worth knowing about for scientific computing, for example the strategy or factory patterns. This style of object-oriented programming is not the one true way to write good code but it's still worth knowing about.

$\endgroup$
4
  • $\begingroup$ Thanks for the references. The Fortran/OOP wikibook is extremely helpful. I'm going to check out the book by Rouson et al. as well. Looks like it'd help me quite a bit. $\endgroup$ – user31765 Nov 15 '20 at 14:40
  • $\begingroup$ I generally recommend reading up on design patterns to anyone willing to listen, but I think that the book you point out (by the "Gang of Four") is a bit outdated. Some of the patterns are just not so relevant any more, in part because they have found their way into standard libraries. $\endgroup$ – Wolfgang Bangerth Nov 16 '20 at 5:02
  • $\begingroup$ @WolfgangBangerth do you know of a better reference? I agree with you that it's pretty old but I haven't found anything better. $\endgroup$ – Daniel Shapero Nov 16 '20 at 17:08
  • $\begingroup$ @DanielShapero I don't, but there are a number of books on design patterns that have appeared over the past few years that are probably worth looking through. I would start with things Kevlin Henney has written. $\endgroup$ – Wolfgang Bangerth Nov 17 '20 at 0:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.