Computation

The main part of the computation involves the implementation of the discretized model by using various Tolosa-lib features. At each time step, several operations need to be managed : the time advancement, the actual model computation, the MPI communications and the timers. The time loop should look like :

do while ( .not. end_time_loop )

    !================================================================================================================!
    !  Time Advancement
    !================================================================================================================!

    [...]

    !================================================================================================================!
    !  Computation
    !================================================================================================================!

    [...]

    !================================================================================================================!
    !  MPI Communications
    !================================================================================================================!

    [...]

end do

Time advancement¶

Firstly, one needs to handle the time advancement. The time is defined by several variables :

nt the time iteration index
tc the current simulation time
dt the time step
ts the total simulation time
end_time_loop is a boolean describing if the total simulation time is reached or not, therefore if the time loop is over or not

Tip

These global variables and other useful tools are defined in the m_common.f90 module.

One has to advance time by adapting the time step if adapt_dt = 1 in the input file, and advance the nt and tc variables. A subroutine advancing nt and tc already exists in the m_common.f90 module.

For example, one can create a subroutine advance_time to manage the time advancement :

SUBROUTINE advance_time

integer(ip) :: ic
real(rp)    :: mes

!================================================================================================================!
!  Adaptive Time Step
!================================================================================================================!

if ( adapt_dt == 1 ) then

    dt = 0._rp

    do ic = 1,mesh%nc

        mes = 2._rp * mesh%cell( ic )%invsurf * mesh%cell( ic )%peri

        dt = max( dt , mes * ( ( .norm. dof%var( ic )%u ) + sqrt( g * dof%var( ic )%h ) ) )

    end do

    dt = cfl / dt

    call mpi_world%min( dt )

end if

!================================================================================================================!
!  Advance nt and tc
!================================================================================================================!

call advance_nt_and_tc( nt , tc , dt )

END SUBROUTINE advance_time

Computation¶

When the time is advanced, one can implement the previously defined Shallow-Water model discretization.

At boundaries¶

Firstly, the computation in the domain's boundaries is handled. One can create a fill_bc subroutine :

SUBROUTINE fill_bc

    integer(ip)   :: ie , iL , iR
    type(primvar) :: sL , sR

    do ie = 1,mesh%neb

        if ( mesh%edgeb( ie )%perio ) cycle

        iL = mesh%edgeb( ie )%cell(1)
        iR = mesh%edgeb( ie )%cell(2)

        sL = dof%var( iL ) .rot. mesh%edgeb( ie )%normal

        select case( mesh%edgeb( ie )%typlim )

            case default

                sR%h   =   sL%h
                sR%u%x = - sL%u%x
                sR%u%y =   sL%u%y

        end select

        dof%var( iR ) = sR .invrot. mesh%edgeb( ie )%normal

    end do

END SUBROUTINE fill_bc

Help on the boundary definition

See Mesh to get all the variables and parameters concerning the boundary conditions to be able to access them.

In the whole domain¶

When the boundary computation is over, one can implement the Euler temporal scheme to be able to compute the Rusanov flux. See Equations to get the schemes' details.

Here is an example of the implementation of these schemes :

SUBROUTINE euler_step_sw

    integer(ip)   :: ie , ic , iL , iR
    type(primvar) :: lflux , nflux , tflux( mesh%nc+mesh%ncdbp ) , sL , sR

    call timer%tic( 'Euler: tflux' )

    tflux(:) = 0._rp

    call timer%tic( 'Euler: loop edge' )

    do ie = 1,mesh%ne

        iL = mesh%edge( ie )%cell(1)
        iR = mesh%edge( ie )%cell(2)

        sL = dof%var( iL ) .rot. mesh%edge( ie )%normal
        sR = dof%var( iR ) .rot. mesh%edge( ie )%normal

        nflux = Rusanov_Flux( sL , sR )

        lflux = ( nflux .invrot. mesh%edge( ie )%normal ) * mesh%edge( ie )%length

                                           tflux( iL ) = tflux( iL ) - lflux
        if ( .not. mesh%edge( ie )%perio ) tflux( iR ) = tflux( iR ) + lflux

    end do

    call timer%tic( 'Euler: loop cell' )

    do ic = 1,mesh%nc

        dof%var( ic )%u = dof%var( ic )%u * dof%var( ic )%h

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

        dof%var( ic )%u = dof%var( ic )%u / dof%var( ic )%h

    end do

END SUBROUTINE euler_step_sw

!**********************************************************************************************************************!
!
!  Rusanov Flux for the SW system of PDE
!
!**********************************************************************************************************************!

ELEMENTAL FUNCTION Rusanov_Flux( sL , sR ) RESULT( flux )

    type(primvar), intent(in) :: sL , sR

    type(primvar) :: flux , fL , fR
    real(rp)      :: lambda

    lambda = max( abs( sL%u%x ) + sqrt( g * sL%h ) , &
             abs( sR%u%x ) + sqrt( g * sR%h ) )

    fL%h   = sL%h * sL%u%x
    fR%h   = sR%h * sR%u%x

    fL%u%x = fL%h * sL%u%x + 0.5_rp * g * sL%h**2
    fL%u%y = fL%h * sL%u%y

    fR%u%x = fR%h * sR%u%x + 0.5_rp * g * sR%h**2
    fR%u%y = fR%h * sR%u%y

    flux%h   = 0.5_rp * ( ( fL%h   + fR%h   ) + lambda * ( sL%h          - sR%h          ) )
    flux%u%x = 0.5_rp * ( ( fL%u%x + fR%u%x ) + lambda * ( sL%h * sL%u%x - sR%h * sR%u%x ) )
    flux%u%y = 0.5_rp * ( ( fL%u%y + fR%u%y ) + lambda * ( sL%h * sL%u%y - sR%h * sR%u%y ) )

END FUNCTION Rusanov_Flux

Info

Some timers have been used in this example, see the Timers subsection to get a general use of them.

MPI communications¶

At the end of the computation, each process has to communicate its primitive variables' values on its subdomain to the other processes to enable them to pursue their computation.

The previously defined unk structure encapsulating the primitive variables values has a subroutine enabling the communication of these variables to the other processes :

call dof%mpicom( mesh )

Timers¶

One can check the time spent on each operation by triggering a timer. In this example, we will check the time spent on each previously described operation. To avoid the definition and initialization of new timers, one can use the tic-toc methods of the timers feature. An example of timers' implementation is the following one :

! Time advancement operations

call timer%tic( 'Advance Time' )
[...]

! Boundary conditions management

call timer%tic( 'Boundary Condition' )
[...]

! Euler time step computation

call timer%tic( 'Euler SW' )
[...]

! MPI communications

call timer%tic( 'MPI Send/Recv' )
[...]

call timer%toc