GET_CFL4CASES.f90

Go to the documentation of this file.

!    Copyright (C) 2012 The SPEED FOUNDATION
!    Author: Ilario Mazzieri
!
!    This file is part of SPEED.
!
!    SPEED is free software; you can redistribute it and/or modify it
!    under the terms of the GNU Affero General Public License as
!    published by the Free Software Foundation, either version 3 of the
!    License, or (at your option) any later version.
!
!    SPEED is distributed in the hope that it will be useful, but
!    WITHOUT ANY WARRANTY; without even the implied warranty of
!    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
!    Affero General Public License for more details.
!
!    You should have received a copy of the GNU Affero General Public License
!    along with SPEED.  If not, see <http://www.gnu.org/licenses/>.
 
!                     'not' deltat could be changed if is extremely small for the time integration
 
      subroutine get_cfl4cases(time_step, nn_loc, loc_n_num, nm, tm, pm, sdeg, &
                             xx_loc, yy_loc, zz_loc, cs_nnz_loc, cs_loc, & 
                             time_step_cfl, fmax, deltat_fixed, &
                             mpi_comm, mpi_np, mpi_id, &
                             ncase,vcase, tcase, zs_elev, zs_all, vs_tria, thick_tria, sub_tag_all, b_failCFL, &
                             damping_type, QS, QP)
          
        use speed_exit_codes
        use speed_par, only: label_testcase
 
        implicit none
 
        include 'SPEED.MPI'
 
        character*3 :: deltat_fixed
 
        integer*4 :: nm, im, nn_loc, cs_nnz_loc, mpi_comm, mpi_np, mpi_err, mpi_id
        integer*4 :: ie, i, j, k, nn, mcode, smcode, sdeg_deltat_cfl, sdeg_npoints
        integer*4 :: damping_type, ncase, icase
        integer*4 :: nel_loc, ic1, ic2, n1, n2, istart, ifin
 
        integer*4, dimension(1) :: pos, pos2
        integer*4, dimension(nm) :: tm, sdeg
        integer*4, dimension(nn_loc) :: loc_n_num
        integer*4, dimension(nn_loc) :: sub_tag_all        
        integer*4, dimension(0:cs_nnz_loc) :: cs_loc
        integer*4, dimension(ncase) ::  tcase, vcase
 
        real*8 :: time_step, time_step_cfl, fmax
        real*8 :: length_min,length,vs_length_min,vs_length,length_vp_min,length_vp
        real*8 :: percent_deltat
        real*8 :: num_of_points,vs_npoints,vp_deltat_cfl
        real*8 :: rho,lambda,mu,vs,vp
        real*8 :: x1_length_vp_min,y1_length_vp_min,z1_length_vp_min
        real*8 :: x2_length_vp_min,y2_length_vp_min,z2_length_vp_min
 
        real*8, dimension(:), allocatable :: ct,ww
        real*8, dimension(:), allocatable :: time_step_glo, vs_length_glo
        real*8, dimension(nn_loc) :: xx_loc, yy_loc, zz_loc
        real*8, dimension(nn_loc) :: zs_elev
        real*8, dimension(nn_loc) :: zs_all
        real*8, dimension(nn_loc) :: vs_tria, thick_tria
        
 
        real*8, dimension(nm,4) :: pm
        real*8, dimension(:,:), allocatable :: dd
        real*8, dimension(nm) :: qs, qp
 
        real*8, dimension(:,:,:), allocatable :: rho_el, lambda_el, mu_el, gamma_el 
 
        logical :: b_failCFL, b_CFL_failure = .false.
 
        nel_loc = cs_loc(0) - 1
 
        length_vp_min = 1.d30
        vs_length_min = 1.d30
        vs_npoints = 1.d30;
        
          
         do ie = 1, nel_loc
 
                im = cs_loc(cs_loc(ie -1) +0)           
                nn = sdeg(im) +1
 
                allocate(rho_el(nn,nn,nn), lambda_el(nn,nn,nn),mu_el(nn,nn,nn),gamma_el(nn,nn,nn)) 
        
                rho_el = pm(im,1) 
                lambda_el = pm(im,2)
                mu_el = pm(im,3)
                gamma_el = pm(im,4)        
            
                do icase = 1, ncase
                 
                  if (vcase(icase) .eq. tm(im)) then
                     call make_eltensor_for_cases(tcase(icase), vcase(icase),&
                                                  nn, rho_el, lambda_el, mu_el, gamma_el,&  
                                                  nn_loc, zs_elev, zs_all, vs_tria, thick_tria,&
                                                  cs_nnz_loc, cs_loc, ie,&
                                                  sub_tag_all, zz_loc, mpi_id, loc_n_num, &
                                                  damping_type, qs(im), qp(im), &
                                                  xx_loc, yy_loc,0,label_testcase)
                  endif
                enddo    
      
 
                smcode=sdeg(im) 
                rho = 0.d0
                lambda = 0.d0
                mu = 0.d0
                      
                do i = 1, nn
                  do j = 1, nn
                     do k = 1, nn
                        rho = rho + rho_el(i,j,k)   
                        lambda = lambda + lambda_el(i,j,k) 
                        mu = mu + mu_el(i,j,k)
                     enddo
                  enddo
                enddo
               
                rho = rho/(nn**3)
                lambda = lambda/(nn**3)
                mu = mu/(nn**3)
                   
                vs = dsqrt(mu/rho)
                vp = dsqrt((lambda+2*mu)/rho)
 
                n1 = nn*nn*(1 -1) +nn*(1 -1) +1
                n2 = nn*nn*(nn -1) +nn*(nn -1) +nn
                istart = cs_loc(ie -1) +n1
                ifin = cs_loc(ie -1) +n2
 
               
               do ic1 = istart, ifin 
                  do ic2 = ic1 + 1, ifin               
               
                  if(cs_loc(ic1) .ne. 0 .and. cs_loc(ic2) .ne. 0) then
               
                     length=dsqrt((xx_loc(cs_loc(ic1))-xx_loc(cs_loc(ic2)))**2 + &
                                  (yy_loc(cs_loc(ic1))-yy_loc(cs_loc(ic2)))**2 + &
                                  (zz_loc(cs_loc(ic1))-zz_loc(cs_loc(ic2)))**2)
                        
                     length_vp = length/vp
 
                    if (length_vp_min .gt. length_vp) then
                        length_vp_min = length_vp
                        x1_length_vp_min = xx_loc(cs_loc(ic1))
                        y1_length_vp_min = yy_loc(cs_loc(ic1))
                        z1_length_vp_min = zz_loc(cs_loc(ic1))
                        x2_length_vp_min = xx_loc(cs_loc(ic2))
                        y2_length_vp_min = yy_loc(cs_loc(ic2))
                        z2_length_vp_min = zz_loc(cs_loc(ic2))
                        vp_deltat_cfl = vp
                        sdeg_deltat_cfl = smcode
                        length_min = length
                     endif
 
 
                     vs_length = vs/length
 
                     if (vs_length_min.gt.vs_length) then
                        vs_length_min = vs_length
                        vs_npoints = vs
                        sdeg_npoints = smcode
                     endif
                                                    
                endif
 
               enddo
             enddo   
      
         deallocate(lambda_el, gamma_el, mu_el, rho_el)
         
        enddo
        
        nn = sdeg_deltat_cfl + 1
 
        time_step_cfl = length_vp_min !/((nn-1)**2)
 
 
        !write(*,*) mpi_id, 'id',  vs_npoints
        !call MPI_BARRIER(mpi_comm, mpi_err)
 
 
        allocate(time_step_glo(mpi_np),vs_length_glo(mpi_np))
        
        call mpi_barrier(mpi_comm, mpi_err)
        call mpi_allgather(time_step_cfl, 1, speed_double, time_step_glo, 1, speed_double, mpi_comm, mpi_err)
        call mpi_allgather(vs_length_min, 1, speed_double, vs_length_glo, 1, speed_double, mpi_comm, mpi_err)
    
        pos  = minloc(time_step_glo)
        pos2 = minloc(vs_length_glo)
 
        !write(*,*) mpi_id, 'id',  pos2(1), vs_length_min, 'oooo',vs_length_glo
        !call MPI_BARRIER(mpi_comm, mpi_err)
         
        deallocate(time_step_glo, vs_length_glo)
 
        if((mpi_id + 1) .eq. pos(1)) then
 
                write(*,'(A)')' '
                write(*,'(A)') '--------Stability for time advancing algorithm --------'
                write(*,'(A,E12.4)') 'Mininum distance beteween grid nodes :', length_vp_min*vp_deltat_cfl
                write(*,'(A,E12.4)') 'Vp for CFL condition         :', vp_deltat_cfl
                write(*,'(A,E12.4,E12.4,E12.4)')'Min. node 1     : ',x1_length_vp_min,y1_length_vp_min,z1_length_vp_min
                write(*,'(A,E12.4,E12.4,E12.4)')'Min. node 2     : ',x2_length_vp_min,y2_length_vp_min,z2_length_vp_min
                write(*,'(A)') '-------------------------------------------------------'
                
                if (time_step.le.time_step_cfl) then
                        percent_deltat=time_step/time_step_cfl*100
                        if (percent_deltat.le.1) then
                                write(*,'(A,E12.4)')'WARNING!!! deltat CFL =',&
                                time_step_cfl
                                write(*,'(A,F6.2,A)')'You are using ',percent_deltat,&
                                        '% of critical time step.'
                                write(*,'(A)')'This time step is excedingly lower the deltat CFL'
                                if (deltat_fixed.eq.'not') then  
                                        write(*,'(A)')'deltat chosen will be substituted with 10% of deltat CFL'
                                        time_step=time_step_cfl*0.1
                                        write(*,'(A,E12.4)')'deltat chosen :',time_step
                                endif
                        elseif (percent_deltat.le.25) then
                                write(*,'(A,E12.4)')'OK!!! deltat CFL =',&
                                time_step_cfl
                                write(*,'(A,F6.2,A)')'You are using ',percent_deltat,&
                                        '% of critical time step.'
                        else
                                write(*,'(A,E12.4)')'WARNING!!! deltat CFL =',&
                                time_step_cfl
                                write(*,'(A,F6.2,A)')'You are using ',percent_deltat,&
                                        '% of critical time step.'
                                write(*,'(A)')'This could be not enough for the correct working of ABC'
                                if (deltat_fixed.eq.'not') then
                                        write(*,'(A)')'deltat chosen will be substituted with 20% of deltat CFL'
                                        time_step=time_step_cfl*0.2
                                        write(*,'(A,E12.4)')'deltat chosen :',time_step
                                endif
                        endif
 
                elseif (time_step .gt. time_step_cfl) then
                        write(*,'(2X,A,E12.4)')'ERROR!!! deltat CFL = ',&
                        time_step_cfl,' < deltat = ',time_step
                        write(*,'(A)')'The time advancing scheme will be unstable!'
                        if (deltat_fixed.eq.'not') then
                                write(*,'(A)')'deltat chosen will be substituted with 20% of deltat CFL'
                                time_step=time_step_cfl*0.2
                                write(*,'(A,E12.4)')'deltat chosen :',time_step
                        else
                                b_cfl_failure = .true.
                        endif
                endif
           
                write(*,'(A)')'-------------------------------------------------------'
                write(*,'(A)')' '
        
                if (b_failcfl .and. b_cfl_failure) then
                        write(*,*) 'CFL does NOT hold, asked to quit. Program finished.'
                        call mpi_finalize(mpi_err)
                        call exit(exit_cfl)
                endif
 
                 !Writing of the results for the maximum number of points per wave length
        
        endif      
 
        call mpi_barrier(mpi_comm, mpi_err)
 
        allocate(ct(nn),ww(nn),dd(nn,nn))
        call make_lgl_nw(nn,ct,ww,dd)
 
        num_of_points = vs_length_min*(ct(1)-ct(nn))/(ct(int(nn/2))-ct(int(nn/2)+1))/fmax
 
        !write(*,*) mpi_id, 'id',  fmax, 'ppppp', vs_length_min, num_of_points
        !call MPI_BARRIER(mpi_comm, mpi_err)
  
        if((mpi_id + 1) .eq. pos2(1)) then
                write(*,'(A)')'-----------Number of points per wave length------------'
                !write(*,'(A,E12.4)')'Max. el. length       :',1/vs_length_min*vs_npoints !length
                write(*,'(A,E12.4)')'Minum vs wrt max frequency  :',vs_npoints
                write(*,'(A,E12.4)')'Points per wavelength       :',num_of_points
                write(*,'(A)') ' '
        endif
 
        deallocate(ct,ww,dd)
        
        call mpi_barrier(mpi_comm, mpi_err)
        
        !stop
        return
               
        end subroutine get_cfl4cases