Assigning pointers to allocatable derived type components in Fortran-CodePudding

I am new to Fortran and am developing a code in which I frequently do operations on two similar derived types (in this case representing models of particles in a scientific computing code). The two types are identical except one has only translational position and velocity components whereas the other also includes angular positions and velocities.

type(particle_A)
  real :: position(3)
  real :: velocity(3)
end type

type(particle_B)
  real :: position(6)
  real :: velocity(6)
end type

I have another derived type which groups arrays of both kinds of particles (together with some other data in the actual application):

type(particleGroup)
  type(particle_A), allocatable :: particleA_array(:)
  integer :: NparticlesA
  type(particle_B), allocatable :: particleB_array(:)
  integer :: NparticlesB
  character(len=30) :: particle_kind
end type

Where both arrays are allocated and re-allocated during runtime if more particles enter the simulation on this particular process. Depending on what particle_kind is set to (during initialization) I want to loop over one of these two arrays and call some subroutine (e.g. for updating the position and velocity). Something like this (in case of looping over the first array)

do iParticle = 1, nParticlesA
  call updateVelocity( particleGroup%particleA_array(iParticle) )
end do

and this for looping over the second array

do iParticle = 1, nParticlesB
  call updateVelocity( particleGroup%particleB_array(iParticle) )
end do

with the generic procedure updateVelocity implemented using an interface block. Is there a way in Fortran to store which array to loop over during initialization (as a pointer?), so that I do not have to check particle_kind at every loop iteration? Or is there some better way of handling this problem?

CodePudding user response：

If you are only interested in using either ParticleA or ParticleB at runtime then you can achieve what you want using a polymorphic array.

You would first have to define a parent class, from which both ParticleA and ParticleB can inherit, something like

module Particle_module
  implicit none

  type, abstract :: Particle
  contains
    procedure(Particle_update_velocity), deferred :: update_velocity
  end type

  interface
    subroutine Particle_update_velocity(this)
      import Particle
      class(Particle), intent(inout) :: this
    end subroutine
  end interface
end module

Then you could write ParticleA as

module ParticleA_module
  use Particle_module
  implicit none

  type, extends(Particle) :: ParticleA
    real :: position(3)
    real :: velocity(3)
  contains
    procedure :: update_velocity => ParticleA_update_velocity
  end type
contains
  subroutine ParticleA_update_velocity(this)
    class(ParticleA), intent(inout) :: this

    ! ParticleA's update_velocity code goes here.
  end subroutine
end module

and analogously for ParticleB.

This would allow you to write ParticleGroup as having a single polymorphic array, which could be either of type ParticleA or ParticleB at runtime, something like

module ParticleGroup_module
  use Particle_module
  implicit none

  type :: ParticleGroup
    class(Particle), allocatable :: particles(:)
    integer :: no_particles
  contains
   procedure :: update_velocities => ParticleGroup_update_velocities
  end type
contains
  subroutine ParticleGroup_update_velocities(this)
    class(ParticleGroup), intent(inout) :: this

    integer :: i

    do i=1,this%no_particles
      call this%particles(i)%update_velocity()
    enddo
  end subroutine
end module

You could then use these classes like

program p
  use Particle_module
  use ParticleA_module
  use ParticleB_module
  use ParticleGroup_module
  implicit none

  type(ParticleGroup) :: particle_group

  ! Fill out the variables of `particle_group`.
  ! Note that the `%particles` array now has the concrete type of `ParticleA`.
  particle_group%particles = [ &
     & ParticleA([0., 0., 0.], [1., 2., 3.]), &
     & ParticleA([1., 1., 1.], [2., 3., 1.]) &
     & ]
  particle_group%no_particles = 2

  ! Call all the `update_velocity` subroutines.
  call particle_group%update_velocities()
end program

In simple examples, keeping particle_group%no_particles separately to size(particle_group%particles) is redundant, but I guess you might need it if you start doing more advanced memory-storage strategies (e.g. C std::vector-style storage).