adds LTS to solver

setup/constants.h.in

@@ -525,12 +525,17 @@
 ! adds overlap elements
   logical,parameter :: LTS_OVERLAP_REGION   = .false.
 
-! uses only 1 single p-level (for debugging)
+! uses only 1 single p-level
   logical,parameter :: LTS_SINGLE_P_LEVEL = .false.
-! uses only 2 p-levels (for debugging)
+! uses only 2 p-levels
   logical,parameter :: LTS_TWO_P_LEVEL = .false.
 
 ! decreases final coarsest time step for stability
+! note: LTS must decrease time step for stability depending on p depth;
+!       to avoid this decrease, one can try to overlap the fine region by one/two elements.
+!       by default, we add an LTS_SAFETY_MARGIN to the suggested LTS time steps and LTS element binning;
+!       where the coarsest time step is calculated as a multiple of the user DT input in Par_file.
+!       here, we can further decrease this coarsest step for LTS simulations to add additional stability.
   logical,parameter :: LTS_DECREASE_DT = .false.
 ! time step margin to use for decrease
   double precision,parameter :: LTS_STABILITY_MARGIN_DT = 0.6d0

src/decompose_mesh/lts_helper.F90

@@ -37,38 +37,18 @@ subroutine lts_add_overlap_elements(nglob,nspec,iglob_p_refine,ispec_p_refine)
   integer,dimension(nspec),intent(inout) :: ispec_p_refine
 
   ! local parameters
-  integer :: ispec,iglob,i,p
+  integer :: ispec,iglob,i,p,icycle
+  ! number of cycles to add overlap by one element layer
+  integer, parameter :: NUMBER_OF_OVERLAP_CYCLES = 1
 
   ! user output
   print *,'  adding overlap region'
+  print *,'    number of overlap cycles = ', NUMBER_OF_OVERLAP_CYCLES
   print *
 
   ! adds overlap by one element
-  ! 1. overlap
-  do ispec = 1,NSPEC
-    ! sets refinement on all points
-    ! this enlarges the finer region by one element
-    p = ispec_p_refine(ispec)
-    do i = 1,NGNOD
-      iglob = elmnts(i,ispec)
-      if ( p > iglob_p_refine(iglob) ) iglob_p_refine(iglob) = p
-    enddo
-  enddo
-  ! re-sets p refinement for element due to overlap addition above
-  ! (adds also elements touching global points with higher refinement)
-  ispec_p_refine(:) = 0
-  do ispec = 1,NSPEC
-    do i = 1,NGNOD
-      iglob = elmnts(i,ispec)
-      p = iglob_p_refine(iglob)
-      if ( p > ispec_p_refine(ispec)) ispec_p_refine(ispec) = p
-    enddo
-  enddo
-  ! checks that all elements have a valid p-value assigned
-  if ( minval(ispec_p_refine(:)) == 0 ) stop 'error ispec_p_refine has zero entry'
-
-  ! 2. overlap
-  if ( .false. ) then
+  do icycle = 1,NUMBER_OF_OVERLAP_CYCLES
+    ! flags all shared nodes according to the highest element p-level
     do ispec = 1,NSPEC
       ! sets refinement on all points
       ! this enlarges the finer region by one element
@@ -90,7 +70,7 @@ subroutine lts_add_overlap_elements(nglob,nspec,iglob_p_refine,ispec_p_refine)
     enddo
     ! checks that all elements have a valid p-value assigned
     if ( minval(ispec_p_refine(:)) == 0 ) stop 'error ispec_p_refine has zero entry'
-  endif
+  enddo
 
   end subroutine lts_add_overlap_elements
 

src/decompose_mesh/lts_setup_elements.F90

@@ -50,8 +50,8 @@ subroutine lts_setup_elements()
 
   ! time steps
   double precision :: dtmin,dtmax
-  double precision :: time_step
   double precision :: deltat_lts
+  double precision :: dt_suggested
 
   real(kind=CUSTOM_REAL),dimension(:),allocatable :: ispec_time_step
 
@@ -155,35 +155,39 @@ subroutine lts_setup_elements()
     distance_min = elemsize * percent_GLL(NGLLX)
 
     ! estimated time step based on CFL condition
-    time_step = COURANT_SUGGESTED * distance_min / max( vp,vs )
+    dt_suggested = COURANT_SUGGESTED * distance_min / max( vp,vs )
+
+    ! cut at a significant number of digits (2 digits)
+    ! example: 0.0734815 -> lpow = (2 - (-1) = 3 -> 0.0730
+    call get_timestep_limit_significant_digit_dp(dt_suggested)
 
     ! stores value
-    ispec_time_step(ispec) = time_step
+    ispec_time_step(ispec) = dt_suggested
 
     ! determines coarsest time step
-    if (dtmax < time_step) then
-      dtmax = time_step
+    if (dtmax < dt_suggested) then
+      dtmax = dt_suggested
     endif
     ! determines finest time step
-    if (dtmin > time_step) then
-      dtmin = time_step
+    if (dtmin > dt_suggested) then
+      dtmin = dt_suggested
     endif
   enddo
-  print *,'  estimated time step min   = ',sngl(dtmin),' seconds'
-  print *,'  estimated time step max   = ',sngl(dtmax),' seconds'
-  print *,'  estimated time step ratio = ',sngl(dtmax/dtmin)
+  print *,'  estimated time step min   = ',dtmin,' seconds'
+  print *,'  estimated time step max   = ',dtmax,' seconds'
+  print *,'  estimated time step ratio = ',dtmax/dtmin
   print *
 
   ! users sets DT = .. in Par_file, we use this as a threshold limit for the minimum time step size
   ! (useful in case CFL is underestimating size)
   if (dtmin > DT) then
-    print *,'  DT time step size set by Par_file: DT = ',sngl(DT),' limits coarsest time step size'
+    print *,'  DT time step size set by Par_file: DT = ',DT,' limits coarsest time step size'
     if (dtmin / DT >= 2.0) then
-      print *,'  Please set DT to higher value    :      ',sngl(dtmin),' otherwise LTS will not be optimal'
+      print *,'  Please set DT to higher value    :      ',dtmin,' otherwise LTS will not be optimal'
     endif
     print *
   endif
-  print *,'  suggested minimum DT time step = ',sngl(dtmin)
+  print *,'  suggested minimum DT time step = ',dtmin
 
   ! we will use now the smallest time step estimate and put elements
   ! into bins with increasing time step size (power of 2 of smallest time step)
@@ -193,17 +197,17 @@ subroutine lts_setup_elements()
   do ispec = 1,nspec
 
     ! gets estimated time step based on CFL condition
-    time_step = ispec_time_step(ispec)
+    dt_suggested = ispec_time_step(ispec)
 
     ! uses a 5% safety margin for binning
-    time_step = ( 1.d0 - LTS_SAFETY_MARGIN ) * time_step
+    dt_suggested = ( 1.d0 - LTS_SAFETY_MARGIN ) * dt_suggested
 
     ! local time step level compared to global time step size
-    p = floor( time_step / dtmin )
+    p = floor( dt_suggested / dtmin )
     if (p < 1) p = 1
 
     ! debug
-    !print *,'ispec',ispec,'time step=',time_step,dtmin,'     p refinement = ',p
+    !print *,'ispec',ispec,'time step=',dt_suggested,dtmin,'     p refinement = ',p
 
     ! debug
     ! uniform distribution for testing

src/decompose_mesh/part_decompose_mesh_hdf5.F90

@@ -1042,7 +1042,7 @@ subroutine write_moho_surface_database_h5(iproc, nspec, &
     endif
   enddo
 
-  !#TODO: check if this works... creates a single group per process
+  !#TODO: HDF5 moho check if this works... creates a single group per process
 
   ! create group name for the database of the current proc
   write(gname_proc,"('proc',i6.6,'_')") iproc

src/decompose_mesh/write_mesh_databases_hdf5.F90

@@ -240,7 +240,7 @@ subroutine write_mesh_databases_hdf5()
 
   if (COUPLE_WITH_INJECTION_TECHNIQUE .or. MESH_A_CHUNK_OF_THE_EARTH) close(124)
 
-  ! #TODO: hdf5 support for fault simulation not implemented yet
+  ! #TODO: HDF5 support for fault simulation not implemented yet
   ! writes fault database
   if (ANY_FAULT) then
     do ipart = 0, nparts-1

src/generate_databases/lts_generate_databases.F90

@@ -98,7 +98,7 @@ subroutine lts_generate_databases()
   ! time steps
   double precision :: dtmin,dtmax
   double precision :: dtmin_glob,dtmax_glob
-  double precision :: time_step
+  double precision :: dt_suggested
 
   real(kind=CUSTOM_REAL),dimension(:),allocatable :: ispec_time_step
 
@@ -215,16 +215,20 @@ subroutine lts_generate_databases()
     !                            NSPEC_AB,NGLOB_AB,ibool,xstore,ystore,zstore)
 
     ! estimated time step based on CFL condition
-    time_step = COURANT_SUGGESTED * distance_min / max( vpmax,vsmax )
+    dt_suggested = COURANT_SUGGESTED * distance_min / max( vpmax,vsmax )
+
+    ! cut at a significant number of digits (2 digits)
+    ! example: 0.0734815 -> lpow = (2 - (-1) = 3 -> 0.0730
+    call get_timestep_limit_significant_digit_dp(dt_suggested)
 
     ! stores value
-    ispec_time_step(ispec) = time_step
+    ispec_time_step(ispec) = dt_suggested
 
     ! determines coarsest time step
-    if (dtmax < time_step) dtmax = time_step
+    if (dtmax < dt_suggested) dtmax = dt_suggested
 
     ! determines finest time step
-    if (dtmin > time_step) dtmin = time_step
+    if (dtmin > dt_suggested) dtmin = dt_suggested
   enddo
 
   ! gets min/max values over all processes
@@ -258,7 +262,7 @@ subroutine lts_generate_databases()
 
   ! user output
   if (myrank == 0) then
-    write(IMAIN,*) '     used minimum DT time step = ',dtmin
+    write(IMAIN,*) '     suggested minimum DT time step = ',dtmin
     call flush_IMAIN()
   endif
 
@@ -270,22 +274,22 @@ subroutine lts_generate_databases()
   do ispec = 1,NSPEC_AB
 
     ! gets estimated time step based on CFL condition
-    time_step = ispec_time_step(ispec)
+    dt_suggested = ispec_time_step(ispec)
 
     ! uses a user-defined safety margin for binning
-    time_step = ( 1.d0 - LTS_SAFETY_MARGIN ) * time_step
+    dt_suggested = ( 1.d0 - LTS_SAFETY_MARGIN ) * dt_suggested
 
     ! absorbing conditions need slightly smaller time steps
     ! if (ABSORBING_CONDITIONS) then
     !   ! uses a 60% safety margin for binning
-    !   time_step = ( 1.d0 - 0.3d0 ) * time_step
+    !   dt_suggested = ( 1.d0 - 0.3d0 ) * dt_suggested
     ! endif
 
     ! local time step level compared to global time step size
-    p = floor( time_step / dtmin )
+    p = floor( dt_suggested / dtmin )
     if (p < 1) p = 1
     ! debug
-    !print *,'ispec',ispec,'time step=',time_step,dtmin,'     p refinement = ',p
+    !print *,'ispec',ispec,'time step=',dt_suggested,dtmin,'     p refinement = ',p
 
     ! debug
     ! uniform distribution for testing
@@ -352,14 +356,14 @@ subroutine lts_generate_databases()
 
   ! user output
   if (myrank == 0) then
-    write(IMAIN,*) '     final lts time step       = ',deltat_lts
+    write(IMAIN,*) '     suggested global coarsest time step       = ',deltat_lts
     write(IMAIN,*)
     call flush_IMAIN()
   endif
   call synchronize_all()
 
   ! re-orders p' -> p = 1 / p'
-  ! note: 1 == coarsest, p_max == smallest such that local time step == dt / p
+  ! note: 1 == coarsest, p_max == finest level such that local time step == dt / p
   do iglob = 1,NGLOB_AB
     p = iglob_p_refine(iglob)
     iglob_p_refine(iglob) = p_max / p
@@ -1251,47 +1255,21 @@ subroutine lts_add_overlap_elements_databases(iglob_p_refine,ispec_p_refine)
   integer,dimension(NSPEC_AB),intent(inout) :: ispec_p_refine
 
   ! local parameters
-  integer :: ispec,iglob,i,j,k,p
+  integer :: ispec,iglob,i,j,k,p,icycle
+  ! number of cycles to add overlap by one element layer
+  integer, parameter :: NUMBER_OF_OVERLAP_CYCLES = 1
 
   ! user output
   if (myrank == 0) then
     write(IMAIN,*) '     adding overlap region'
+    write(IMAIN,*) '       number of overlap cycles = ', NUMBER_OF_OVERLAP_CYCLES
     write(IMAIN,*)
     call flush_IMAIN()
   endif
 
   ! adds overlap by one element
-  ! 1. overlap
-  do ispec = 1,NSPEC_AB
-    ! sets refinement on all points
-    ! this enlarges the finer region by one element
-    p = ispec_p_refine(ispec)
-    do k = 1,NGLLZ
-      do j = 1,NGLLY
-        do i = 1,NGLLX
-          iglob = ibool(i,j,k,ispec)
-          if (p > iglob_p_refine(iglob)) iglob_p_refine(iglob) = p
-        enddo
-      enddo
-    enddo
-  enddo
-  ! re-sets p refinement for element due to overlap addition above
-  ! (adds also elements touching global points with higher refinement)
-  ispec_p_refine(:) = 0
-  do ispec = 1,NSPEC_AB
-    do k = 1,NGLLZ
-      do j = 1,NGLLY
-        do i = 1,NGLLX
-          iglob = ibool(i,j,k,ispec)
-          p = iglob_p_refine(iglob)
-          if (p > ispec_p_refine(ispec)) ispec_p_refine(ispec) = p
-        enddo
-      enddo
-    enddo
-  enddo
-
-  ! 2. overlap
-  if (.false.) then
+  do icycle = 1,NUMBER_OF_OVERLAP_CYCLES
+    ! flags all shared nodes according to the highest element p-level
     do ispec = 1,NSPEC_AB
       ! sets refinement on all points
       ! this enlarges the finer region by one element
@@ -1319,7 +1297,7 @@ subroutine lts_add_overlap_elements_databases(iglob_p_refine,ispec_p_refine)
         enddo
       enddo
     enddo
-  endif
+  enddo
 
   end subroutine lts_add_overlap_elements_databases
 

src/inverse_problem_for_model/rules.mk

@@ -163,6 +163,10 @@ inverse_problem_for_model_OBJECTS += \
 	$O/locate_point.spec.o \
 	$O/locate_receivers.spec.o \
 	$O/locate_source.spec.o \
+	$O/lts_assemble_MPI_vector.spec.o \
+	$O/lts_global_step.spec.o \
+	$O/lts_iterate_time.spec.o \
+	$O/lts_newmark_update.spec.o \
 	$O/lts_setup.spec.o \
 	$O/make_gravity.spec.o \
 	$O/noise_tomography.spec.o \

src/inverse_problem_for_model/specfem_interface/specfem_interface_mod.F90

@@ -1570,7 +1570,12 @@ subroutine  CheckModelSuitabilityForModeling(ModelIsSuitable)
      ! suggested timestep
      vel_max = max( vpmax,vsmax )
      dt_suggested = COURANT_SUGGESTED * distance_min / vel_max
-     dt_suggested_glob = min( dt_suggested_glob, dt_suggested)
+
+     ! cut at a significant number of digits (2 digits)
+     ! example: 0.0734815 -> lpow = (2 - (-1) = 3 -> 0.0730
+     call get_timestep_limit_significant_digit(dt_suggested)
+
+     dt_suggested_glob = min(dt_suggested_glob, dt_suggested)
 
   enddo
 

src/shared/check_mesh_resolution.f90

@@ -294,7 +294,13 @@ subroutine check_mesh_resolution(NSPEC_AB,NGLOB_AB,ibool,xstore,ystore,zstore, &
 
     ! suggested timestep
     dt_suggested = COURANT_SUGGESTED * distance_min / vel_max
-    dt_suggested_glob = min( dt_suggested_glob, dt_suggested)
+
+    ! cut at a significant number of digits (2 digits)
+    ! example: 0.0734815 -> lpow = (2 - (-1) = 3 -> 0.0730
+    call get_timestep_limit_significant_digit(dt_suggested)
+
+    ! global suggested timestep
+    dt_suggested_glob = min(dt_suggested_glob, dt_suggested)
 
     ! debug: for vtk output
     if (SAVE_MESH_FILES) tmp2(ispec) = pmax
@@ -736,7 +742,7 @@ end subroutine check_mesh_resolution
 !
 !    ! suggested timestep
 !    dt_suggested = COURANT_SUGGESTED * distance_min / max( vpmax,vp2max,vsmax )
-!    dt_suggested_glob = min( dt_suggested_glob, dt_suggested)
+!    dt_suggested_glob = min(dt_suggested_glob, dt_suggested)
 !
 !    ! debug: for vtk output
 !    if (SAVE_MESH_FILES) tmp2(ispec) = pmax