SPEED.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/>.
 
 
! Here starts the code SPEED
 
      program speed
 
      use max_var
      use str_mesh 
      use str_mesh_scratch                   
      use dgjump
      use speed_par
      use speed_par_dg
 
      implicit none
 
      include 'SPEED.MPI'
      
 
!*****************************************************************************************************************
      
      allocate(mpi_stat(speed_status_size))
      mpi_comm = speed_comm
      
      call initialization(mpi_comm,mpi_np,mpi_id,mpi_ierr)
 
      start = mpi_wtime()
      
      if (mpi_id.eq.0) then
         write(*,'(A)')''
         write(*,'(A)')''
         write(*,'(A)')'*******************************************************'
         write(*,'(A)')'*                                                     *'
         write(*,'(A)')'*                        SPEED                        *'
         write(*,'(A)')'*          SPectral Elements in Elastodynamics        *'
         write(*,'(A)')'*              with Discontinuous Galerkin            *'
         write(*,'(A)')'*                                                     *'
         write(*,'(A)')©'*           PoliMi, 2012, All Rights Reserved        *'
         write(*,'(A)')'*                                                     *'
         write(*,'(A)')'*******************************************************'
         write(*,'(A)')''
         write(*,'(A)')''
      endif
 
!***********************************************************************************
!     READ INPUT FILES AND ALLOCATE VARIABLES
!***********************************************************************************
 
 
     call read_input_files()
 
     
!***********************************************************************************
!     PARTITION OF THE GRID AND GENERATION OF LOCAL CONNECTIVITY 
!***********************************************************************************
 
 
     allocate(elem_domain(nelem));  elem_domain = mpi_id;
     call make_partition_and_mpi_files()
 
 
!***********************************************************************************
!    START SETUP 4 MPI COMUNICATION (ONLY BETWEEN CONFORMING REGIONS)
!***********************************************************************************
 
 
     call make_setup_mpi_conforming()
 
                   
!***********************************************************************************
!    NEW SPECTRAL LOCAL CONNECTIVITY
!***********************************************************************************
 
 
     call make_spx_grid_with_loc_numeration()
 
 
!*****************************************************************************************
!    BUILD SEISMIC MOMENT
!*****************************************************************************************
 
      if (srcmodflag.eq.0) then
          call make_seismic_moment_or_explosive_source()
      elseif (srcmodflag.eq.1) then
          call make_seismic_moment_pointsources_nhfault()
      endif
 
 
!*****************************************************************************************
!    DAMPING 
!*****************************************************************************************
      
     make_damping_yes_or_not = 0
      
     if (damping_type .eq. 1) then !Kosloff & Kosloff damping
        do im = 1,nmat        
           if (abs(prop_mat(im,4)) .gt. 10e-5) then 
              make_damping_yes_or_not = 1                
           endif                                        
        enddo                                
     
        if(mpi_id.eq.0 .and. make_damping_yes_or_not .eq. 1) then        
                write(*,'(A)')                         
                write(*,'(A)') '----------------Building the Damping matrix---------------'
        endif
        
        if(mpi_id.eq.0 .and. make_damping_yes_or_not .eq. 0) then        
                write(*,'(A)')                                                        
                write(*,'(A)')'----------------Building the Damping matrix---------------'                        
                write(*,'(A)')'ATT: There are no materials with damping defined on'        
        endif
        
     elseif(damping_type .eq. 2) then !Standard Linear Solid damping
        if(mpi_id.eq.0) then   
                write(*,'(A)')                     
                write(*,'(A)') '------------Computing anelastic coefficients-----------'
        endif
        !Assume 3 SLS for a range ~ [fmax/10, fmax*10] Hz
        allocate(y_lambda(nmat,n_sls),y_mu(nmat,n_sls),frequency_range(n_sls))
        
        y_lambda=0.d0;  y_mu=0.d0;         
                 
        if (fmax .ne. 0.d0) &                
            call make_anelastic_coefficients(nmat, n_sls, prop_mat, qs, qp, fmax, &
                                         y_lambda, y_mu, frequency_range, mpi_id)
           
            ! Check coefficient signs
            if (b_failoncoeffs) then
              do im = 1,  nmat
                 do i = 1, n_sls
                   if (y_lambda(im,i) .lt. 0.d0 .or. y_mu(im,i) .lt. 0) then
                      call mpi_finalize(mpi_ierr)
                      if (mpi_id .eq. 0) write(*,'(A)') 'NEGATIVE anelastic coefficient found, program ended.'
                      if (mpi_ierr.ne.0) write(*,'(A,I6)')'MPI Finalization error - proc : ',mpi_id
                      call exit(exit_anelastic)
                   endif
                 enddo
 
              enddo
 
            endif
     elseif(damping_type .eq. 3)  then  !Rayleigh damping
        if(mpi_id.eq.0) then   
                write(*,'(A)')                     
                write(*,'(A)') '------------Computing Rayleigh coefficients-----------'
        endif
        allocate(a0_ray(nmat),a1_ray(nmat))
        a0_ray = 1.d0; a1_ray = 1.d0;
        
              
        if (fmax .ne. 0.d0) &                
            call make_rayleigh_coefficients(nmat, a0_ray, a1_ray, prop_mat, qs, fmax,mpi_id)
         
     endif   
      
 
!*****************************************************************************************
!   Not Honoring Enhanced (Reading Material properties from Tomo File)          
!*****************************************************************************************
    if (nmat_nhe .gt. 0) then
        call make_nh_enhanced()
        call mpi_barrier(mpi_comm, mpi_ierr)
    endif
 
 
!***********************************************************************************
!     FIND OSCILLATOR POSITION
!***********************************************************************************
   if (sdofflag.eq.1) then
      call read_system_position()
   endif
!*****************************************************************************************
!   BUILD LOAD MATRIX AND SEISMIC MOMENT TENSOR OR EXPLOSIVE SOURCE
!*****************************************************************************************
 
 
   call make_load_matrix()
 
      
!**************************************************************************************            
!   FIND MONITORED NODES  (PGM & MLST input files) + deallocate(elem_domain)
!**************************************************************************************
 
 
     call find_monitor_position()
     deallocate(elem_domain);   !Should we deallocate the elem_domain?
 
 
!*****************************************************************************************
!   CASES 4 NOT HONORING          
!*****************************************************************************************
 
   
    if(n_case .gt. 0)  then
                                      
       if (mpi_id.eq.0) write(*,'(A)')'----------------Making not-honoring case---------------'                  
 
       allocate(zs_elev(nnod_loc), zs_all(nnod_loc), sub_tag_all(nnod_loc), vs_tria(nnod_loc), thick(nnod_loc))  
       
       if(mpi_id .eq.0) then
          start1=mpi_wtime()  
       endif 
       
       call make_nothonoring(local_node_num, nnod_loc,&
                             n_case, tag_case, val_case, tol_case, &
                             con_nnz_loc, con_spx_loc, nmat, tag_mat, sdeg_mat, &                        
                             xx_spx_loc, yy_spx_loc, zz_spx_loc, zs_elev, zs_all, vs_tria, thick, &
                             sub_tag_all, mpi_id)
 
  
 
                                                                            
      ncase = 1
      CALL mpi_barrier(mpi_comm, mpi_ierr) 
      
      if (mpi_id.eq.0) then
         start2=mpi_wtime()
         write(*,'(A)')
         write(*,'(A,F20.3,A)')'Set-up time CASE = ',start2-start1,' s'
         write(*,'(A)')'Made.'                  
      endif                        
                                    
    else
    
      ncase = 0
 
    endif
 
!*****************************************************************************************
!   RANDOMIZED PARAMETERS          
!*****************************************************************************************
    if (nmat_rnd .gt. 0) then
        allocate(lambda_rnd(nnod_loc),mu_rnd(nnod_loc),rho_rnd(nnod_loc))
        
        call make_random_param(local_node_num, nnod_loc, &
                               nmat_rnd, rand_mat,  &
                               xx_spx_loc, yy_spx_loc, zz_spx_loc, &
                               lambda_rnd, mu_rnd, rho_rnd, mpi_id)
 
    endif
 
 
!*****************************************************************************************
!    SETUP FOR BOUNDARY CONDITIONS    
!*****************************************************************************************
 
     
     call make_boundary_conditions()
 
 
!*****************************************************************************************
!    SETUP FOR DG - INTERFACE CONDITIONS          
!*****************************************************************************************
 
     
     call make_dg_interface_conditions()
 
 
!***********************************************************************************************      
!    COMPUTE CFL CONDITION
!***********************************************************************************************
 
      if (mpi_id.eq.0) write(*,'(A)') 
      if (mpi_id.eq.0) write(*,'(A)')'------------Setting calculation parameters-------------'
 
 
      
      if(n_case .gt. 0) then 
 
        
          call get_cfl4cases(deltat, nnod_loc, local_node_num, &
                      nmat, tag_mat, prop_mat, sdeg_mat,&
                      xx_spx_loc, yy_spx_loc, zz_spx_loc, &
                      con_nnz_loc, con_spx_loc, &
                      deltat_cfl, fmax, deltat_fixed, mpi_comm, mpi_np, mpi_id, &
                      n_case,val_case, tag_case, zs_elev, zs_all, vs_tria, thick, sub_tag_all, b_failcfl, &
                      damping_type, qs, qp)
 
      else if (nmat_nhe.gt.0) then
          call get_cfl4nhe(deltat, nnod_loc, local_node_num, &
                             nmat, tag_mat, prop_mat, sdeg_mat, &
                             xx_spx_loc, yy_spx_loc, zz_spx_loc, &
                             con_nnz_loc, con_spx_loc, & 
                             rho_nhe, lambda_nhe, mu_nhe, &
                             deltat_cfl, fmax, deltat_fixed, &
                             mpi_comm, mpi_np, mpi_id, b_failcfl)
      else
       
          call get_cfl(deltat, nnod_loc, local_node_num, &
                      nmat, tag_mat, prop_mat, sdeg_mat,&
                      xx_spx_loc, yy_spx_loc, zz_spx_loc, &
                      con_nnz_loc, con_spx_loc, &
                      deltat_cfl, fmax, deltat_fixed, mpi_comm, mpi_np, mpi_id, b_failcfl) 
       endif
 
 
      nts = int((tstop -tstart) / deltat)
      tstop = tstart + dfloat(nts)*deltat
      if (mpi_id.eq.0) then
         write(*,'(A,I15)')  'Number of time-steps : ',nts
         write(*,'(A,E14.5)')'Start time           : ',tstart
         write(*,'(A,E14.5)')'Final time           : ',tstop
      endif
      
!**************************************************************************************      
!     SETUP FOR SNAPSHOTS
!**************************************************************************************      
      
!      if (nsnaps.gt.0) then
!        do i = 1, nsnaps
!          itersnap(i) = int((tsnap(i) - tstart) / deltat)
!          if (itersnap(i) .lt. 0) itersnap(i) = 0
!          if (itersnap(i) .gt. nts) itersnap(i) = nts
!        enddo
        
!        initial_snap = 1  
!        if(tstart .gt. 0.d0) then
!           allocate(set_initial_snap(nsnaps))
!           set_initial_snap = tstart - tsnap 
           
!           do i = 2, nsnaps
!                   if(set_initial_snap(i) .lt. set_initial_snap(i-1) .and. set_initial_snap(i) .ge. 0) then
!                      initial_snap = i
!                   endif   
!           enddo        
!           deallocate(set_initial_snap)
!         endif                 
!      endif
 
!***********************************************************************************
!    WRITE MESH VTK
!***********************************************************************************
      !call WRITE_VTK_MESH(nnod_loc, local_node_num, nmat, tag_mat, prop_mat, &
      !sdeg_mat, xx_spx_loc, yy_spx_loc, zz_spx_loc, con_nnz_loc, con_spx_loc, mpi_id)
 
!**************************************************************************************      
!    END SETUP 
!**************************************************************************************      
 
      finish = mpi_wtime()
      time_in_seconds = finish-start
      
      if (mpi_id.eq.0) then
         write(*,'(A)') 
         write(*,'(A)')'-------------------------------------------------------'
         write(*,'(A,F20.4,A)')'Set-up time = ',time_in_seconds,' s'
         write(*,'(A)')'-------------------------------------------------------'
         write(*,'(A)')
      endif
 
      if (b_setuponly) then
        call mpi_finalize(mpi_ierr)
        if (mpi_ierr.ne.0) write(*,'(A,I6)')'MPI Finalization error - proc : ',mpi_id
        if (mpi_id .eq. 0) then
            write(*, '(A)') 'Set-up ended. Program finished.'
        endif
        call exit(exit_setup)
      endif
 
!**************************************************************************************      
!    TIME INTEGRATION  
!**************************************************************************************      
 
      if (mpi_id.eq.0) then
         write(*,'(A)')
         write(*,'(A)')'------------Beginning of the time-loop-----------------'
         write(*,'(A)')
      endif
           
     
     call time_loop(el_new)
 
!**************************************************************************************      
!     FINALIZE AND DEALLOCATE MEMORY
!**************************************************************************************      
 
      call mpi_finalize(mpi_ierr)
      
      if (mpi_ierr.ne.0) write(*,'(A,I6)')'MPI Finalization error - proc : ',mpi_id
                     
      if ((mpi_id.eq.0).and.(opt_out_data.gt.0)) then
         write(*,'(A)');  write(*,'(A)') 'Print output'      
      endif
            
      call deallocate_variables()
 
      if (mpi_id.eq.0) write(*,'(A)')
      if (mpi_id.eq.0) write(*,'(A)')'Bye.'
      
      call exit(exit_normal)
            
      end program speed