Computation Loop

This page describes the main time-stepping loop for the Shallow-Water equations.

Main Loop Structure

! Time loop
t = 0.0_rp
iter = 0

do while (t < t_end)

  ! 1. Compute time step
  call compute_dt(dt, cfl)

  ! 2. Time integration (RK scheme)
  call time_step_rk(h, u, v, dt)

  ! 3. Apply boundary conditions
  call apply_bc(h, u, v)

  ! 4. Update time
  t = t + dt
  iter = iter + 1

  ! 5. Output if needed
  if (mod(iter, output_freq) == 0) then
    call write_output(t, h, u, v)
  end if

  ! 6. Display progress
  if (rank == 0 .and. mod(iter, 100) == 0) then
    print '(A,F8.2,A,I8)', 'Time: ', t, 's  Iteration: ', iter
  end if

end do

Time Step Computation

The time step must satisfy the CFL condition:

subroutine compute_dt(dt, cfl)
  real(rp), intent(out) :: dt
  real(rp), intent(in) :: cfl
  real(rp) :: c, dx_min
  integer :: ic

  dt = 1.0e10_rp

  ! Loop over cells
  do ic = 1, mesh%nc
    ! Cell size
    dx_min = mesh%cell_size(ic)

    ! Wave speed
    c = sqrt(u(ic)**2 + v(ic)**2) + sqrt(g * h(ic))

    ! Local time step
    dt = min(dt, cfl * dx_min / c)
  end do

  ! MPI reduction for minimum
  call mpi_allreduce_min(dt)

end subroutine compute_dt

Time Integration - Runge-Kutta

RK2 Scheme

subroutine time_step_rk2(h, u, v, dt)
  real(rp), intent(inout) :: h(:), u(:), v(:)
  real(rp), intent(in) :: dt
  real(rp), allocatable :: h1(:), u1(:), v1(:)

  ! Allocate intermediate arrays
  allocate(h1(mesh%nc), u1(mesh%nc), v1(mesh%nc))

  ! Stage 1
  call compute_rhs(h, u, v, rhs_h, rhs_u, rhs_v)
  h1 = h + dt * rhs_h
  u1 = u + dt * rhs_u
  v1 = v + dt * rhs_v

  ! Stage 2
  call compute_rhs(h1, u1, v1, rhs_h, rhs_u, rhs_v)
  h = 0.5_rp * (h + h1 + dt * rhs_h)
  u = 0.5_rp * (u + u1 + dt * rhs_u)
  v = 0.5_rp * (v + v1 + dt * rhs_v)

  deallocate(h1, u1, v1)

end subroutine time_step_rk2

RK4 Scheme

For higher accuracy:

subroutine time_step_rk4(h, u, v, dt)
  ! Four-stage Runge-Kutta
  ! Implementation similar to RK2 with 4 stages
end subroutine time_step_rk4

Right-Hand Side Computation

Compute spatial derivatives and fluxes:

subroutine compute_rhs(h, u, v, rhs_h, rhs_u, rhs_v)
  real(rp), intent(in) :: h(:), u(:), v(:)
  real(rp), intent(out) :: rhs_h(:), rhs_u(:), rhs_v(:)
  integer :: ic, ie, ic1, ic2
  real(rp) :: flux_h, flux_u, flux_v

  ! Initialize
  rhs_h = 0.0_rp
  rhs_u = 0.0_rp
  rhs_v = 0.0_rp

  ! Start timer
  call t_flux%start()

  ! Loop over edges
  do ie = 1, mesh%ne
    ic1 = mesh%edge_cells(1, ie)
    ic2 = mesh%edge_cells(2, ie)

    ! Compute numerical flux
    call compute_edge_flux(ie, h, u, v, flux_h, flux_u, flux_v)

    ! Accumulate to cells
    rhs_h(ic1) = rhs_h(ic1) - flux_h / mesh%cell_vol(ic1)
    rhs_h(ic2) = rhs_h(ic2) + flux_h / mesh%cell_vol(ic2)

    rhs_u(ic1) = rhs_u(ic1) - flux_u / mesh%cell_vol(ic1)
    rhs_u(ic2) = rhs_u(ic2) + flux_u / mesh%cell_vol(ic2)

    rhs_v(ic1) = rhs_v(ic1) - flux_v / mesh%cell_vol(ic1)
    rhs_v(ic2) = rhs_v(ic2) + flux_v / mesh%cell_vol(ic2)
  end do

  ! Stop timer
  call t_flux%stop()

  ! Add source terms
  call add_source_terms(h, u, v, rhs_h, rhs_u, rhs_v)

end subroutine compute_rhs

Numerical Flux Computation

HLL Flux

subroutine compute_hll_flux(ie, h, u, v, flux_h, flux_u, flux_v)
  integer, intent(in) :: ie
  real(rp), intent(in) :: h(:), u(:), v(:)
  real(rp), intent(out) :: flux_h, flux_u, flux_v

  integer :: ic1, ic2
  real(rp) :: h1, h2, u1, u2, v1, v2
  real(rp) :: un1, un2, c1, c2
  real(rp) :: sl, sr, sm
  real(rp) :: nx, ny, length

  ! Get cells
  ic1 = mesh%edge_cells(1, ie)
  ic2 = mesh%edge_cells(2, ie)

  ! Get edge normal
  nx = mesh%edge_normal(1, ie)
  ny = mesh%edge_normal(2, ie)
  length = mesh%edge_length(ie)

  ! Left and right states
  h1 = h(ic1); h2 = h(ic2)
  u1 = u(ic1); u2 = u(ic2)
  v1 = v(ic1); v2 = v(ic2)

  ! Normal velocities
  un1 = u1*nx + v1*ny
  un2 = u2*nx + v2*ny

  ! Wave speeds
  c1 = sqrt(g * h1)
  c2 = sqrt(g * h2)

  sl = min(un1 - c1, un2 - c2)
  sr = max(un1 + c1, un2 + c2)

  ! HLL flux
  if (sl >= 0.0_rp) then
    flux_h = h1 * un1
    flux_u = h1 * u1 * un1 + 0.5_rp * g * h1**2 * nx
    flux_v = h1 * v1 * un1 + 0.5_rp * g * h1**2 * ny
  else if (sr <= 0.0_rp) then
    flux_h = h2 * un2
    flux_u = h2 * u2 * un2 + 0.5_rp * g * h2**2 * nx
    flux_v = h2 * v2 * un2 + 0.5_rp * g * h2**2 * ny
  else
    ! Mix left and right states
    sm = (sr * flux_r - sl * flux_l) / (sr - sl)
    ! ... (complete HLL formulation)
  end if

  ! Multiply by edge length
  flux_h = flux_h * length
  flux_u = flux_u * length
  flux_v = flux_v * length

end subroutine compute_hll_flux

Source Terms

Bottom Friction

subroutine add_friction(h, u, v, rhs_u, rhs_v)
  real(rp), intent(in) :: h(:), u(:), v(:)
  real(rp), intent(inout) :: rhs_u(:), rhs_v(:)
  integer :: ic
  real(rp) :: umag, cf

  do ic = 1, mesh%nc
    umag = sqrt(u(ic)**2 + v(ic)**2)
    cf = g * manning**2 / h(ic)**(4.0_rp/3.0_rp)

    rhs_u(ic) = rhs_u(ic) - cf * umag * u(ic)
    rhs_v(ic) = rhs_v(ic) - cf * umag * v(ic)
  end do

end subroutine add_friction

Boundary Conditions

subroutine apply_bc(h, u, v)
  real(rp), intent(inout) :: h(:), u(:), v(:)
  integer :: ib, ic

  do ib = 1, mesh%nb
    ic = mesh%boundary_cells(ib)

    select case (mesh%bc_type(ib))
      case (BC_WALL)
        ! Reflect normal velocity
        call apply_wall_bc(ib, u(ic), v(ic))

      case (BC_OPEN)
        ! Radiation condition
        call apply_open_bc(ib, h(ic), u(ic), v(ic))

      case (BC_INFLOW)
        ! Prescribed values
        call apply_inflow_bc(ib, h(ic), u(ic), v(ic))
    end select
  end do

end subroutine apply_bc

Performance Optimization

MPI Communication

! Exchange ghost cell data
call mesh%exchange_data(h)
call mesh%exchange_data(u)
call mesh%exchange_data(v)

OpenMP Parallelization

!$OMP PARALLEL DO
do ic = 1, mesh%nc
  ! Compute cell contribution
end do
!$OMP END PARALLEL DO

Next Steps

• Finalization - Output and cleanup
• Tolosa-sw - Full production model