Initialization

The following paragraphs will suggest an initialization of Tolosa-lib concerning the Shallow water equations. A general Tolosa-lib initialization will be presented, as well as several implementations of structures of primitive variables.

Tolosa initialization

To initialize Tolosa-lib, one has to call the Init_Tolosa routine. Here, the mesh and input are specified to initialize the mesh of the domain as well as the input parameters written in the input.txt file.

type(msh) :: mesh
type(cli) :: input
[...]

call Init_Tolosa( input=input , mesh=mesh )

call input%add_get( 'g' , default = '9.81' , kind=rtype , value=g )

Simply by calling this routine, the MPI environment, the arguments passed to the command line or specified in the input.txt file, the timers, and the mesh are initialized. Here the gravity constant can be defined in the input.txt file, therefore, this argument has been added to the ones to be read.

Primitive variables

To be able to get the primitive variables from the previously detailed Shallow Water discretization, one can choose to either create an array of structures (AoS) of primitive variables or a structure of arrays (SoA) of primitive variables. The two methods will be detailed in the following paragraphs.

Either way, the primitive variables will be stored in an unk structure and allocated and initialized with :

[...]
type(unk) :: dof

[...]

call dof%init( mesh )

Array of structures (AoS)

AoS is the more conventional layout, in which data for different fields is interleaved. This is often more intuitive, and supported directly by most programming languages.¹

Source

The source for this example of AoS use can be found in the ./Tolosa-lib/examples/sw/ Tolosa-lib directory (See Download to get the Tolosa-lib sources).

TYPE primvar

  real(rp) :: h
  type(vec2d) :: u

CONTAINS

  procedure, pass( self ) :: equal_primvar
  procedure, pass( self ) :: equal_primvar_scalar
  procedure               :: add_primvar
  procedure               :: sub_primvar
  procedure, pass( self ) :: scalar_primvar
  procedure, pass( self ) :: primvar_scalar
  procedure               :: primvar_div_scalar
  procedure, pass( self ) :: rot_primvar
  procedure, pass( self ) :: invrot_primvar

  generic :: assignment(=) => equal_primvar , &
                            equal_primvar_scalar

  generic :: operator(+) => add_primvar
  generic :: operator(-) => sub_primvar
  generic :: operator(*) => scalar_primvar , &
                          primvar_scalar
  generic :: operator(/) => primvar_div_scalar

  generic :: operator(.rot.)    => rot_primvar
  generic :: operator(.invrot.) => invrot_primvar

END TYPE primvar

Here, a structure has been created to get the primitive variables and . This structure stores the height value in h and the speed values in u a vec2d object (See 2D vectors). This structure gathers these primitive variables’ values at a specific point of the domain. Therefore, an AoS of this structure is created to get the primitive variables’ values at each cell of the domain.

!===================================================================================================================!
!  Array of Structures of Primitive Variables (AoS)
!===================================================================================================================!

TYPE unk

  type(primvar), allocatable :: var(:)

CONTAINS

  procedure, pass( self ) ::    init =>    init_dof
  procedure, pass( self ) ::   alloc =>   alloc_dof
  procedure, pass( self ) :: dealloc => dealloc_dof
  procedure, pass( self ) ::  mpicom =>     com_dof

END TYPE unk

Here, var(:) is an array of primvar structures. The dimension of this array will be the number of mesh cells to enable a structure of primitive variables to be linked to each mesh cell.

One can easily understand that the implementation of the Shallow Water discretization will be specific to the use of AoS (see the detailed implementation in the source code).

type(msh) :: mesh
type(unk) :: dof

[...]

call dof%init( mesh )

[...]

  do ic = 1,mesh%nc

    [...]
    dof%var( ic ) = dof%var( ic ) + tflux( ic ) * dt * mesh%cell( ic )%invsurf
    [...]

  end do

[...]

call dof%dealloc

Structure of Arrays (SoA)

SoA is a layout separating elements of a record into one parallel array per field. The manipulation is easier with packed SIMD instructions in most instruction set architectures. If only a specific part of the record is needed, only those parts need to be iterated over, allowing more data to fit onto a single cache line.¹

Source

The source for this example of SoA use can be found in the ./Tolosa-lib/examples/sw_soa/ Tolosa-lib directory (See Download to get the Tolosa-lib sources).

!===================================================================================================================!
!  Structure of Primitive Variables Arrays (SoA)
!===================================================================================================================!

TYPE unk

  real(rp), allocatable :: h(:)
  real(rp), allocatable :: u(:)
  real(rp), allocatable :: v(:)

CONTAINS

  procedure, pass( dof ) ::    init =>    init_dof
  procedure, pass( dof ) ::   alloc =>   alloc_dof
  procedure, pass( dof ) :: dealloc => dealloc_dof
  procedure, pass( dof ) ::  mpicom =>     com_dof

END TYPE unk

Unlike the previous method, this one contains a set of arrays for each primitive variable. This unk structure is a Structure of Arrays (SoA). Each array encapsulates the concerned variable distibution on the domain (at each mesh cell).

One can easily understand that the implementation of the Shallow Water discretization will be specific to the use of SoA (see the detailed implementation in the source code).

type(msh) :: mesh
type(unk) :: dof

[...]

call dof%init( mesh )

[...]

  do ic = 1,mesh%nc

    [...]
    dof%h( ic ) = dof%h( ic ) + tflux_h( ic ) * dt * mesh%cell( ic )%invsurf
    dof%u( ic ) = dof%u( ic ) + tflux_u( ic ) * dt * mesh%cell( ic )%invsurf
    dof%v( ic ) = dof%v( ic ) + tflux_v( ic ) * dt * mesh%cell( ic )%invsurf
    [...]

  end do

[...]

call dof%dealloc

Note

Here, the computation of the Rusanov flux still uses the primvar type since the implementation is more intuitive. However, the fluxes’ records are saved in separate arrays.

SoA and mesh structure copy

This method also uses SoA, and instead of accessing the whole mesh structure, it copies elements of the structure in local arrays.

[...]

integer(ip), allocatable :: edge_cell(:,:)
type(vec2d), allocatable :: normal(:)
real(rp)   , allocatable :: invsurf(:) , peri(:) , length(:)
logical(lp), allocatable :: perio(:)

[...]

allocate( edge_cell( 2 , mesh%ne ) )
allocate( normal   (     mesh%ne ) )
allocate( length   (     mesh%ne ) )
allocate( perio    (     mesh%ne ) )

do ie = 1,mesh%ne

  edge_cell(1,ie) = mesh%edge( ie )%cell(1)
  edge_cell(2,ie) = mesh%edge( ie )%cell(2)
  normal   (  ie) = mesh%edge( ie )%normal
  length   (  ie) = mesh%edge( ie )%length
  perio    (  ie) = mesh%edge( ie )%perio

end do

allocate( invsurf( mesh%nc ) )
allocate( peri   ( mesh%nc ) )

do ic = 1,mesh%nc

  invsurf( ic ) = mesh%cell( ic )%invsurf
  peri   ( ic ) = mesh%cell( ic )%peri

end do

[...]

Starting timer(s)

When the initialization of the primitive variables, mesh, inputs, timers, etc. is completed, one can start the main timer by calling :

call timer%start( timerid=timer_main )

Any other timer can be initialized and started if needed.

Footnotes

1. AoS and SoA ↩ ↩²