Merge branch 'master' into reconnection

src/driver.F90

@@ -40,16 +40,13 @@ program amun
   use equations     , only : initialize_equations, finalize_equations
   use evolution     , only : initialize_evolution, finalize_evolution
   use evolution     , only : advance, new_time_step
-  use evolution     , only : n, t, dt
+  use evolution     , only : step, time, dt
   use integrals     , only : init_integrals, clear_integrals, store_integrals
   use interpolations, only : initialize_interpolations, finalize_interpolations
   use io            , only : initialize_io, write_data, write_restart_data     &
-                           , restart_job
+                           , read_restart_data, next_tout
   use mesh          , only : initialize_mesh, finalize_mesh
   use mesh          , only : generate_mesh, store_mesh_stats
-#ifdef MPI
-  use mesh          , only : redistribute_blocks
-#endif /* MPI */
   use mpitools      , only : initialize_mpitools, finalize_mpitools
   use mpitools      , only : setup_mpi
 #ifdef MPI
@@ -82,6 +79,12 @@ program amun
   logical, dimension(3) :: per = .true.
   integer               :: nmax  = 0, ndat = 1, nres = -1, ires = -1
   real                  :: tmax  = 0.0d+00, trun = 9.999d+03, tsav = 3.0d+01
+  real(kind=8)          :: dtnext = 0.0d+00
+
+! flag to adjust time precisely to the snapshots
+!
+  logical        , save :: precise_snapshots = .false.
+  character(len=255)    :: prec_snap         = "off"
 
 ! temporary variables
 !
@@ -255,6 +258,23 @@ program amun
 !
   trun   = trun - tsav / 6.0d+01
 
+! initialize dtnext
+!
+  dtnext = 2.0d+00 * tmax
+
+! get the precise snapshot flag
+!
+  call get_parameter_string ("precise_snapshots", prec_snap)
+
+! set the precise snapshot flag
+!
+  if (prec_snap == "on"  ) precise_snapshots = .true.
+  if (prec_snap == "ON"  ) precise_snapshots = .true.
+  if (prec_snap == "true") precise_snapshots = .true.
+  if (prec_snap == "TRUE") precise_snapshots = .true.
+  if (prec_snap == "yes" ) precise_snapshots = .true.
+  if (prec_snap == "YES" ) precise_snapshots = .true.
+
 ! get integral calculation interval
 !
   call get_parameter_integer("ndat" , ndat)
@@ -363,6 +383,14 @@ program amun
 !
   call initialize_problems()
 
+! initialize boundaries module and print info
+!
+  if (master) then
+    write (*,*)
+    write (*,"(1x,a)"         ) "Snapshots:"
+    write (*,"(4x,a22,1x,'='1x,a)") "precise snapshot times", trim(prec_snap)
+  end if
+
 ! initialize module IO
 !
   call initialize_io()
@@ -375,10 +403,6 @@ program amun
 !
     call initialize_mesh(nrun, master, iret)
 
-! initialize the integrals module
-!
-    call init_integrals(.true.)
-
 ! generate the initial mesh, refine that mesh to the desired level according to
 ! the initialized problem
 !
@@ -386,33 +410,34 @@ program amun
 
 ! store mesh statistics
 !
-    call store_mesh_stats(n, t)
+    call store_mesh_stats(step, time)
 
 ! calculate new timestep
 !
-    call new_time_step()
+    call new_time_step(dtnext)
+
+! initialize the integrals module
+!
+    call init_integrals(.true.)
 
   else
 
+! increase the run number
+!
+    nrun  = nres + 1
+
 ! initialize the mesh module
 !
     call initialize_mesh(nrun, master, iret)
 
+! reconstruct the meta and data block structures from a given restart file
+!
+    call read_restart_data()
+
 ! initialize the integrals module
 !
     call init_integrals(.false.)
 
-! reconstruct the meta and data block structures from a given restart file
-!
-    call restart_job()
-
-#ifdef MPI
-! redistribute blocks between processors in case the number of processors has
-! changed
-!
-    call redistribute_blocks()
-#endif /* MPI */
-
   end if
 
 #ifdef MPI
@@ -460,7 +485,7 @@ program amun
 ! get current time in seconds
 !
   if (master) &
-    tbeg = t
+    tbeg = time
 
 ! print progress info on master processor
 !
@@ -482,32 +507,36 @@ program amun
 #if defined INTEL || defined PATHSCALE
     write(*,'(i8,2(1x,1pe14.6),2x,i8,2x,1i4.1,"d",1i2.2,"h",1i2.2,"m"' //      &
             ',1i2.2,"s",15x,a1,$)')                                            &
-                              n, t, dt, get_nleafs(), ed, eh, em, es, char(13)
+                        step, time, dt, get_nleafs(), ed, eh, em, es, char(13)
 #else /* INTEL | PATHSCALE */
     write(*,'(i8,2(1x,1pe14.6),2x,i8,2x,1i4.1,"d",1i2.2,"h",1i2.2,"m"' //      &
             ',1i2.2,"s",15x,a1)',advance="no")                                 &
-                              n, t, dt, get_nleafs(), ed, eh, em, es, char(13)
+                        step, time, dt, get_nleafs(), ed, eh, em, es, char(13)
 #endif /* INTEL | PATHSCALE */
 
   end if
 
 ! main loop
 !
-  do while((nsteps < nmax) .and. (t <= tmax) .and. (iterm == 0))
+  do while((nsteps <= nmax) .and. (time < tmax) .and. (iterm == 0))
+
+! get the next snapshot time
+!
+    if (precise_snapshots) dtnext = next_tout() - time
 
 ! performe one step evolution
 !
-    call advance()
+    call advance(dtnext)
 
 ! advance the iteration number and time
 !
-    t      = t + dt
+    time   = time   + dt
-    n      = n + 1
+    step   = step   +  1
     nsteps = nsteps +  1
 
 ! store mesh statistics
 !
-    call store_mesh_stats(n, t)
+    call store_mesh_stats(step, time)
 
 ! store integrals
 !
@@ -535,7 +564,7 @@ program amun
 
 ! calculate days, hours, seconds
 !
-      ec   = int(tm_curr * (tmax - t) / max(1.0e-8, t - tbeg), kind = 4)
+      ec   = int(tm_curr * (tmax - time) / max(1.0e-8, time - tbeg), kind = 4)
       es   = max(0, int(mod(ec, 60)))
       em   = int(mod(ec / 60, 60))
       eh   = int(ec / 3600)
@@ -546,11 +575,11 @@ program amun
 #if defined INTEL || defined PATHSCALE
       write(*,'(i8,2(1x,1pe14.6),2x,i8,2x,1i4.1,"d",1i2.2,"h",1i2.2,"m"' //    &
               ',1i2.2,"s",15x,a1,$)')                                          &
-                              n, t, dt, get_nleafs(), ed, eh, em, es, char(13)
+                        step, time, dt, get_nleafs(), ed, eh, em, es, char(13)
 #else /* INTEL | PATHSCALE */
       write(*,'(i8,2(1x,1pe14.6),2x,i8,2x,1i4.1,"d",1i2.2,"h",1i2.2,"m"' //    &
               ',1i2.2,"s",15x,a1)',advance="no")                               &
-                              n, t, dt, get_nleafs(), ed, eh, em, es, char(13)
+                        step, time, dt, get_nleafs(), ed, eh, em, es, char(13)
 #endif /* INTEL | PATHSCALE */
 
     end if

src/evolution.F90

@@ -50,8 +50,8 @@ module evolution
 
 ! time variables
 !
-  integer, save :: n       = 0
+  integer, save :: step    = 0
-  real   , save :: t       = 0.0d+00
+  real   , save :: time    = 0.0d+00
   real   , save :: dt      = 1.0d+00
   real   , save :: dtn     = 1.0d+00
 
@@ -66,7 +66,7 @@ module evolution
 
 ! declare public variables
 !
-  public :: cfl, n, t, dt, dtn
+  public :: cfl, step, time, dt, dtn
 
 !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 !
@@ -201,10 +201,13 @@ module evolution
 !   Subroutine advances the solution by one time step using the selected time
 !   integration method.
 !
+!   Arguments:
+!
+!     dtnext - next time step;
 !
 !===============================================================================
 !
-  subroutine advance()
+  subroutine advance(dtnext)
 
 ! include external procedures
 !
@@ -218,12 +221,16 @@ module evolution
 ! local variables are not implicit by default
 !
     implicit none
+
+! input variables
+!
+    real, intent(in) :: dtnext
 !
 !-------------------------------------------------------------------------------
 !
 ! find new time step
 !
-    call new_time_step()
+    call new_time_step(dtnext)
 
 ! advance the solution using the selected method
 !
@@ -252,6 +259,124 @@ module evolution
   end subroutine advance
 !
 !===============================================================================
+!
+! subroutine NEW_TIME_STEP:
+! ------------------------
+!
+!   Subroutine estimates the new time step from the maximum speed in the system
+!   and source term constraints.
+!
+!   Arguments:
+!
+!     dtnext - next time step;
+!
+!===============================================================================
+!
+  subroutine new_time_step(dtnext)
+
+! include external procedures
+!
+    use equations     , only : maxspeed, cmax, cmax2
+#ifdef MPI
+    use mpitools      , only : reduce_maximum_real, reduce_maximum_integer
+#endif /* MPI */
+
+! include external variables
+!
+    use blocks        , only : block_data, list_data
+    use coordinates   , only : adx, ady, adz
+    use coordinates   , only : toplev
+
+! local variables are not implicit by default
+!
+    implicit none
+
+! input variables
+!
+    real, intent(in) :: dtnext
+
+! local pointers
+!
+    type(block_data), pointer :: pblock
+
+! local variables
+!
+    integer                   :: iret
+    integer(kind=4)           :: lev
+    real                      :: cm, dx_min
+
+! local parameters
+!
+    real, parameter           :: eps = tiny(cmax)
+!
+!-------------------------------------------------------------------------------
+!
+! reset the maximum speed, and the highest level
+!
+    cmax   = eps
+    lev    = 1
+
+! iterate over all data blocks in order to find the maximum speed among them
+! and the highest level which is required to obtain the spatial step
+!
+    pblock => list_data
+    do while (associated(pblock))
+
+! find the maximum level occupied by blocks (can be smaller than toplev)
+!
+      lev = max(lev, pblock%meta%level)
+
+! obtain the maximum speed for the current block
+!
+      cm = maxspeed(pblock%q(:,:,:,:))
+
+! compare global and local maximum speeds
+!
+      cmax = max(cmax, cm)
+
+! assiociate the pointer with the next block
+!
+      pblock => pblock%next
+
+    end do
+
+#ifdef MPI
+! find maximum speed in the system from all processors
+!
+    call reduce_maximum_real   (cmax, iret)
+    call reduce_maximum_integer(lev , iret)
+#endif /* MPI */
+
+! calculate squared cmax
+!
+    cmax2 = cmax * cmax
+
+! find the smallest spatial step
+!
+#if NDIMS == 2
+    dx_min = min(adx(lev), ady(lev))
+#endif /* NDIMS == 2 */
+#if NDIMS == 3
+    dx_min = min(adx(lev), ady(lev), adz(lev))
+#endif /* NDIMS == 3 */
+
+! calcilate the new time step
+!
+    dtn = dx_min / max(cmax, eps)
+
+! calculate the new time step
+!
+    dt  = cfl * dtn
+
+! round the time
+!
+    if (dtnext > 0.0d+00) dt = min(dt, dtnext)
+
+!-------------------------------------------------------------------------------
+!
+  end subroutine new_time_step
+!
+!===============================================================================
 !!
 !!***  PRIVATE SUBROUTINES  ****************************************************
 !!
@@ -599,113 +724,6 @@ module evolution
 !-------------------------------------------------------------------------------
 !
   end subroutine update_variables
-!
-!===============================================================================
-!
-! subroutine NEW_TIME_STEP:
-! ------------------------
-!
-!   Subroutine estimates the new time step from the maximum speed in the system
-!   and source term constraints.
-!
-!
-!===============================================================================
-!
-  subroutine new_time_step()
-
-! include external procedures
-!
-    use equations     , only : maxspeed, cmax, cmax2
-#ifdef MPI
-    use mpitools      , only : reduce_maximum_real, reduce_maximum_integer
-#endif /* MPI */
-
-! include external variables
-!
-    use blocks        , only : block_data, list_data
-    use coordinates   , only : adx, ady, adz
-    use coordinates   , only : toplev
-
-! local variables are not implicit by default
-!
-    implicit none
-
-! local pointers
-!
-    type(block_data), pointer :: pblock
-
-! local variables
-!
-    integer                   :: iret
-    integer(kind=4)           :: lev
-    real                      :: cm, dx_min
-
-! local parameters
-!
-    real, parameter           :: eps = tiny(cmax)
-!
-!-------------------------------------------------------------------------------
-!
-! reset the maximum speed, and the highest level
-!
-    cmax   = eps
-    lev    = 1
-
-! iterate over all data blocks in order to find the maximum speed among them
-! and the highest level which is required to obtain the spatial step
-!
-    pblock => list_data
-    do while (associated(pblock))
-
-! find the maximum level occupied by blocks (can be smaller than toplev)
-!
-      lev = max(lev, pblock%meta%level)
-
-! obtain the maximum speed for the current block
-!
-      cm = maxspeed(pblock%q(:,:,:,:))
-
-! compare global and local maximum speeds
-!
-      cmax = max(cmax, cm)
-
-! assiociate the pointer with the next block
-!
-      pblock => pblock%next
-
-    end do
-
-#ifdef MPI
-! find maximum speed in the system from all processors
-!
-    call reduce_maximum_real   (cmax, iret)
-    call reduce_maximum_integer(lev , iret)
-#endif /* MPI */
-
-! calculate squared cmax
-!
-    cmax2 = cmax * cmax
-
-! find the smallest spatial step
-!
-#if NDIMS == 2
-    dx_min = min(adx(lev), ady(lev))
-#endif /* NDIMS == 2 */
-#if NDIMS == 3
-    dx_min = min(adx(lev), ady(lev), adz(lev))
-#endif /* NDIMS == 3 */
-
-! calcilate the new time step
-!
-    dtn = dx_min / max(cmax, eps)
-
-! calculate the new time step
-!
-    dt  = cfl * dtn
-
-!-------------------------------------------------------------------------------
-!
-  end subroutine new_time_step
 
 !===============================================================================
 !

src/integrals.F90

@@ -138,7 +138,7 @@ module integrals
     use blocks   , only : block_meta, block_data, list_data
     use coordinates, only : ib, ie, jb, je, kb, ke
     use coordinates, only : advol
-    use evolution, only : n, t, dt
+    use evolution, only : step, time, dt
     use mpitools , only : master
 #ifdef MPI
     use mpitools , only : reduce_sum_real_array
@@ -211,7 +211,7 @@ module integrals
 ! close integrals.dat
 !
     if (master) then
-      write(funit,"(i8,12(1x,1pe15.8))") n, t, dt, arr(1:10)
+      write(funit,"(i8,12(1x,1pe15.8))") step, time, dt, arr(1:10)
     end if
 !
 !-------------------------------------------------------------------------------

src/makefile

@@ -141,7 +141,7 @@ integrals.o      : integrals.F90 blocks.o coordinates.o equations.o \
 interpolations.o : interpolations.F90 blocks.o coordinates.o parameters.o      \
                    timers.o
 io.o             : io.F90 blocks.o coordinates.o equations.o error.o           \
-                   evolution.o mpitools.o random.o refinement.o timers.o
+                   evolution.o mesh.o mpitools.o random.o refinement.o timers.o
 mesh.o           : mesh.F90 blocks.o coordinates.o domains.o equations.o       \
                    error.o interpolations.o mpitools.o problems.o refinement.o \
                    timers.o

src/mesh.F90

@@ -1402,7 +1402,7 @@ module mesh
 
 ! local buffer for data block exchange
 !
-    real(kind=8)   , dimension(nv,im,jm,km) :: rbuf
+    real(kind=8)   , dimension(2,nv,im,jm,km) :: rbuf
 #endif /* MPI */
 
 !-------------------------------------------------------------------------------
@@ -1451,7 +1451,8 @@ module mesh
 
 ! copy data to buffer
 !
-            rbuf(:,:,:,:) = pmeta%data%u(:,:,:,:)
+            rbuf(1,:,:,:,:) = pmeta%data%u(:,:,:,:)
+            rbuf(2,:,:,:,:) = pmeta%data%q(:,:,:,:)
 
 ! send data
 !
@@ -1480,7 +1481,8 @@ module mesh
 
 ! coppy the buffer to data block
 !
-            pmeta%data%u(:,:,:,:) = rbuf(:,:,:,:)
+            pmeta%data%u(:,:,:,:) = rbuf(1,:,:,:,:)
+            pmeta%data%q(:,:,:,:) = rbuf(2,:,:,:,:)
 
           end if ! nproc == n