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 \(\begin{pmatrix} h \\ u \\ v \end{pmatrix}\) 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.1
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).
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 | |
Here, a structure has been created to get the primitive variables \(h\) and \(\mathbf{u}\). 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.
32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 | |
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.1
Source
The source for this example of AoS use can be found in the ./Tolosa-lib/examples/sw_soa/ Tolosa-lib directory (See Download to get the Tolosa-lib sources).
32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 | |
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 is 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.