Skip to content
Tolosa Documentation

Computation Loop

This page describes the main time-stepping loop for the Shallow-Water equations example. The source can be found in the ./Tolosa-lib/examples/sw/ directory (See Download).

Structure

Section titled “Structure”

The loop runs until end_time_loop is set to .true. by the time advancement routine. Inside the loop:

do while ( .not. end_time_loop )


   call timer%tic( 'Advance Time' )
   call advance_time


   call timer%tic( 'Boundary Condition' )
   call fill_bc


   call euler_step_sw


   call timer%tic( 'MPI Send/Recv' )
   call dof%mpicom( mesh )
   call timer%toc


end do

Time Advancement

Section titled “Time Advancement”

The advance_time subroutine computes the adaptive time step (when adapt_dt = 1) and advances the time iteration index nt and simulation time tc.

The CFL-based time step is computed as:

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 )

where mes is an approximation of the inverse of the local mesh size, and mpi_world%min reduces the time step across all MPI processes. The final time loop advancement is done by the Tolosa-lib routine advance_nt_and_tc.

Boundary Conditions

Section titled “Boundary Conditions”

The fill_bc subroutine fills the ghost cells at physical boundaries. For each boundary edge (looping on mesh%neb), it rotates the left state into the local frame, applies the BC, and stores the result in the right ghost cell:

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

Euler Time Step

Section titled “Euler Time Step”

The euler_step_sw subroutine computes the Rusanov flux on each edge and applies the forward Euler update. For each interior edge:

  1. Rotate both states onto the local normal frame using .rot.
  2. Compute the Rusanov flux
  3. Rotate the flux back to the global frame using .invrot.
  4. Accumulate the flux contribution (scaled by edge length) into the two neighboring cells
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

The cell update then applies the flux balance divided by cell area:

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

MPI Communication

Section titled “MPI Communication”

After the Euler update, the halo cells of neighboring MPI subdomains are refreshed by calling:

call dof%mpicom( mesh )

which sends and receives the h, u%x, and u%y fields across subdomain boundaries (see Computation Loop source for the mpicom implementation).

Next Steps

Section titled “Next Steps”