Add functionality to control time step in af_advance

The flux routine returns the maximal wave speed, and the advance routine has an
extra argument dt_lim (maximal time step)

examples/Makefile

@@ -10,12 +10,11 @@ include $(AF_DIR)/src/makerules.make
 PROGS_XD := random_refinement poisson_basic poisson_benchmark advection		\
  computational_domain boundary_conditions poisson_neumann particles_to_grid	\
  ghostcell_benchmark particles_gravity poisson_coarse_solver implicit_diffusion	\
- poisson_dielectric poisson_helmholtz reaction_diffusion dielectric_surface
+ poisson_dielectric poisson_helmholtz reaction_diffusion dielectric_surface advection_v2
 
 PROGS_2D := $(PROGS_XD:%=%_2d) poisson_cyl poisson_cyl_dielectric		\
 simple_streamer poisson_cyl_analytic poisson_helmholtz_cyl h_domain_expts h_bc	\
-h_advection solid_body_rotation amr_solid_body_rotation euler_gas_dynamics	\
-advection_v2
+h_advection solid_body_rotation amr_solid_body_rotation euler_gas_dynamics
 
 PROGS_3D := $(PROGS_XD:%=%_3d) poisson_div_cleaning
 

examples/advection_v2.f90

@@ -16,12 +16,12 @@ program advection
 
   type(af_t)         :: tree
   type(ref_info_t)   :: refine_info
-  integer            :: refine_steps, time_steps, output_cnt
-  integer            :: n, n_steps
-  real(dp)           :: dt, time, end_time, err, sum_err2
+  integer            :: refine_steps, output_cnt
+  real(dp)           :: dt, dt_lim, time, time_prev_refine
+  real(dp)           :: end_time, err, sum_err2
   real(dp)           :: sum_phi, sum_phi_t0
-  real(dp)           :: dt_adapt, dt_output
-  real(dp)           :: velocity(NDIM), dr_min(NDIM)
+  real(dp)           :: dt_amr, dt_output
+  real(dp)           :: velocity(NDIM)
   character(len=100) :: fname
 
   print *, "Running advection_" // DIMNAME // ""
@@ -43,11 +43,11 @@ program advection
        [DTIMES(box_size)], &
        periodic=[DTIMES(.true.)])
 
-  output_cnt = 0
-  time       = 0
-  dt_adapt   = 0.01_dp
-  dt_output  = 0.5_dp
-  end_time   = 5.0_dp
+  output_cnt  = 0
+  time        = 0
+  dt_amr      = 0.01_dp
+  dt_output   = 0.5_dp
+  end_time    = 5.0_dp
   velocity(:) = 0.0_dp
   velocity(1) = 1.0_dp
   velocity(2) = -1.0_dp
@@ -78,19 +78,14 @@ program advection
   ! Fill ghost cells for variables i_phi on the sides of all boxes
   call af_gc_tree(tree, [i_phi])
 
-  time_steps = 0
-
   call af_tree_sum_cc(tree, i_phi, sum_phi_t0)
 
+  dt = 0.5_dp / (sum(abs(velocity/af_lvl_dr(tree, tree%highest_lvl))) + &
+       epsilon(1.0_dp))
+  time_prev_refine = time
+
   ! Starting simulation
   do
-     time_steps = time_steps + 1
-     dr_min  = af_min_dr(tree)
-     dt      = 0.5_dp / (sum(abs(velocity/dr_min)) + epsilon(1.0_dp))
-
-     n_steps = ceiling(dt_adapt/dt)
-     dt      = dt_adapt / n_steps
-
      if (output_cnt * dt_output <= time) then
         output_cnt = output_cnt + 1
         write(fname, "(A,I0)") "advection_v2_" // DIMNAME // "_", output_cnt
@@ -115,13 +110,15 @@ program advection
 
      if (time > end_time) exit
 
-     do n = 1, n_steps
-        call af_advance(tree, dt, time, [i_phi], &
-             af_midpoint_method, forward_euler)
-     end do
+     call af_advance(tree, dt, dt_lim, time, [i_phi], &
+          af_midpoint_method, forward_euler)
+     dt = 0.8_dp * dt_lim
 
-     call af_gc_tree(tree, [i_phi])
-     call af_adjust_refinement(tree, refine_routine, refine_info, 1)
+     if (time > time_prev_refine + dt_amr) then
+        call af_gc_tree(tree, [i_phi])
+        call af_adjust_refinement(tree, refine_routine, refine_info, 1)
+        time_prev_refine = time
+     end if
   end do
 
 contains
@@ -210,17 +207,22 @@ contains
     end select
   end function solution
 
-  subroutine forward_euler(tree, dt, time, s_deriv, s_prev, s_out)
+  subroutine forward_euler(tree, dt, dt_lim, time, s_deriv, s_prev, s_out)
     type(af_t), intent(inout) :: tree
     real(dp), intent(in)      :: dt
     real(dp), intent(in)      :: time
+    real(dp), intent(out)     :: dt_lim
     integer, intent(in)       :: s_deriv
     integer, intent(in)       :: s_prev
     integer, intent(in)       :: s_out
+    real(dp)                  :: wmax(NDIM)
 
-    call flux_generic_tree(tree, 1, [i_phi+s_deriv], [i_flux], &
+    call flux_generic_tree(tree, 1, [i_phi+s_deriv], [i_flux], wmax, &
          max_wavespeed, get_flux)
     call flux_update_densities(tree, dt, 1, [i_phi], [i_flux], s_prev, s_out)
+
+    ! Compute maximal time step
+    dt_lim = 1.0_dp / sum(wmax/af_lvl_dr(tree, tree%highest_lvl))
   end subroutine forward_euler
 
   subroutine max_wavespeed(n_values, n_var, flux_dim, u, w)

src/m_af_advance.f90

@@ -15,10 +15,11 @@ module m_af_advance
   integer, parameter :: req_copies(n_integrators) = [1, 2, 2]
 
   interface
-     subroutine subr_feuler(tree, dt, time, s_deriv, s_prev, s_out)
+     subroutine subr_feuler(tree, dt, dt_lim, time, s_deriv, s_prev, s_out)
        import
        type(af_t), intent(inout) :: tree
        real(dp), intent(in)      :: dt      !< Time step
+       real(dp), intent(out)     :: dt_lim  !< Computed time step limit
        real(dp), intent(in)      :: time    !< Current time
        integer, intent(in)       :: s_deriv !< State to compute derivatives from
        integer, intent(in)       :: s_prev  !< Previous state
@@ -31,11 +32,12 @@ module m_af_advance
 contains
 
   !> Compute generic finite volume flux
-  subroutine af_advance(tree, dt, time, i_cc, time_integrator, forward_euler)
+  subroutine af_advance(tree, dt, dt_lim, time, i_cc, time_integrator, forward_euler)
     type(af_t), intent(inout) :: tree
-    real(dp), intent(in)      :: dt
-    real(dp), intent(inout)   :: time
-    integer, intent(in)       :: i_cc(:)
+    real(dp), intent(in)      :: dt      !< Current time step
+    real(dp), intent(out)     :: dt_lim  !< Time step limit
+    real(dp), intent(inout)   :: time    !< Current time
+    integer, intent(in)       :: i_cc(:) !< Index of cell-centered variables
     integer, intent(in)       :: time_integrator
     procedure(subr_feuler)    :: forward_euler
 
@@ -47,16 +49,16 @@ contains
 
     select case (time_integrator)
     case (af_forward_euler)
-       call forward_euler(tree, dt, time, 0, 0, 0)
+       call forward_euler(tree, dt, dt_lim, time, 0, 0, 0)
        time = time + dt
     case (af_midpoint_method)
-       call forward_euler(tree, 0.5_dp * dt, time, 0, 0, 1)
-       call forward_euler(tree, dt, time, 1, 0, 0)
+       call forward_euler(tree, 0.5_dp * dt, dt_lim, time, 0, 0, 1)
+       call forward_euler(tree, dt, dt_lim, time, 1, 0, 0)
        time = time + dt
     case (af_heuns_method)
-       call forward_euler(tree, dt, time, 0, 0, 1)
+       call forward_euler(tree, dt, dt_lim, time, 0, 0, 1)
        time = time + dt
-       call forward_euler(tree, dt, time, 1, 1, 1)
+       call forward_euler(tree, dt, dt_lim, time, 1, 1, 1)
        call combine_substeps(tree, i_cc, 2, [0, 1], [0.5_dp, 0.5_dp], 0)
     end select
 

src/m_af_flux_schemes.f90

@@ -280,7 +280,7 @@ contains
     !$omp end parallel
   end subroutine flux_update_densities
 
-  subroutine flux_generic_tree(tree, n_vars, i_cc, i_flux, &
+  subroutine flux_generic_tree(tree, n_vars, i_cc, i_flux, wmax, &
        max_wavespeed, flux_from_primitives, to_primitive, to_conservative)
     use m_af_restrict
     use m_af_core
@@ -288,6 +288,7 @@ contains
     integer, intent(in)            :: n_vars         !< Number of variables
     integer, intent(in)            :: i_cc(n_vars)   !< Cell-centered variables
     integer, intent(in)            :: i_flux(n_vars) !< Flux variables
+    real(dp), intent(out)          :: wmax(NDIM)     !< Maximum wave speed found
     !> Compute the maximum wave speed
     procedure(subr_max_wavespeed)  :: max_wavespeed
     !> Compute the flux from primitive variables
@@ -302,13 +303,15 @@ contains
     ! Ensure ghost cells near refinement boundaries can be properly filled
     call af_restrict_ref_boundary(tree, i_cc)
 
-    !$omp parallel private(lvl, i)
+    wmax(:) = 0
+
+    !$omp parallel private(lvl, i) reduction(max:wmax)
     do lvl = 1, tree%highest_lvl
        !$omp do
        do i = 1, size(tree%lvls(lvl)%leaves)
           call flux_generic_box(tree, tree%lvls(lvl)%leaves(i), tree%n_cell, &
-               n_vars, i_cc, i_flux, max_wavespeed, flux_from_primitives, &
-               to_primitive, to_conservative)
+               n_vars, i_cc, i_flux, wmax, max_wavespeed, &
+               flux_from_primitives, to_primitive, to_conservative)
        end do
        !$omp end do
     end do
@@ -320,7 +323,7 @@ contains
   end subroutine flux_generic_tree
 
   !> Compute generic finite volume flux
-  subroutine flux_generic_box(tree, id, nc, n_vars, i_cc, i_flux, &
+  subroutine flux_generic_box(tree, id, nc, n_vars, i_cc, i_flux, wmax, &
        max_wavespeed, flux_from_primitives, to_primitive, to_conservative)
     use m_af_types
     use m_af_ghostcell
@@ -330,6 +333,7 @@ contains
     integer, intent(in)            :: n_vars         !< Number of variables
     integer, intent(in)            :: i_cc(n_vars)   !< Cell-centered variables
     integer, intent(in)            :: i_flux(n_vars) !< Flux variables
+    real(dp), intent(inout)        :: wmax(NDIM)     !< Maximum wave speed found
     !> Compute the maximum wave speed
     procedure(subr_max_wavespeed)  :: max_wavespeed
     !> Compute the flux from primitive variables
@@ -397,8 +401,12 @@ contains
                 call to_conservative(nc+1, n_vars, u_r)
              end if
 
+             w_l = max(w_l, w_r) ! Get maximum of left/right wave speed
              call flux_kurganovTadmor_1d(nc+1, n_vars, flux_l, flux_r, &
-                  u_l, u_r, max(w_l, w_r), flux)
+                  u_l, u_r, w_l, flux)
+
+             ! Store maximum wave speed
+             wmax(flux_dim) = max(wmax(flux_dim), maxval(w_l))
 
              ! Store the computed fluxes
              select case (flux_dim)