# Best practice for ADTs in computational science with Fortran

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.