test.f90 Source File

program test

  use parameters
  use time_info
  use DG_structures
  use grid
  use scalar_schw
  use all_integrators
  use observers
  use geodesic_schwarzschild
  use osculating_schwarzschild
  use world_tube
  use scalar_schw_eff
  use self_force_observer
  use self_force_base
  use numerics, only : time_window
  use analytic_circular_orbit
  use singular_observer

  implicit none

  class(equation), pointer  :: eqp
  type(scal_schw), allocatable, target :: my_eq
  integer(ip) :: i, j, k
  type(ref_element) :: relem
  class(equation_pointer), dimension(:), allocatable :: eqs
  type(cobserver), dimension(:,:), allocatable :: cobs
  type(sing_observer) :: sobs
  class(equation), pointer :: orbit
  complex(wp), dimension(:,:), allocatable :: psi, dpsidt, dpsidr
  real(wp) :: r, phi, ur, En, Lz, tfac, dtfac_dt, dt2fac_dt2

  call read_parameters ()

  print*,'p_orb = ', p_orb
  print*,'ecc = ', ecc
  relem = ref_element( order )

  call init_grid_coordinates ( relem )

  call init_time ( t_initial )

  if (use_world_tube) then
    wt = wtube()

!    call wt%win%output(rho)
!    call wt%dwin%output(rho)
!    call wt%d2win%output(rho)
!    call wt%boundary_info%output(rho)
  end if

  if ( use_generic_orbit ) then
    if ( use_osculating_orbit ) then
      allocate ( osc_schw :: orbit )
    else
      allocate ( geod_schw :: orbit )
    end if
  else
    allocate ( circular_orbit :: orbit )
  end if 
  call orbit%init ( )

  print*,'rho%vname = ', rho%vname
  print*,'equation_name = ', equation_name
  print*,'n_elems = ', n_elems
  print*,'order = ', order
  allocate(my_eq)
  call my_eq%init ( )

  if ( output_coords_for_exact ) then
    call my_eq%output_coords ( )
  end if

  call sobs%init ( (/ p_orb /), my_eq%r_star, my_eq%effs )

  allocate(eqs(2))
  eqs(1)%p => my_eq
  eqs(2)%p => orbit
!  allocate(eqs(2))
!  eqs(1)%p => orbit
  print*,'my_eq%eq_data(1,1)%vname = ', my_eq%eq_data(1,1)%vname
  print*,'my_eq%eq_data(2,1)%vname = ', my_eq%eq_data(2,1)%vname
  print*,'my_eq%eq_data(3,1)%vname = ', my_eq%eq_data(3,1)%vname

  if (use_particle) then
    call time_window ( get_current_time (), tsigma, torder, tfac, &
                                     dtfac_dt, dt2fac_dt2 )
    print*,'Setting time window: ', get_current_time ( ), tfac, dtfac_dt, dt2fac_dt2
! Tell the effective source about the time 
    if ( turn_on_source_smoothly ) then
      call my_eq%effs%set_time_window ( tfac, dtfac_dt, dt2fac_dt2, -1 )
    else
      call my_eq%effs%set_time_window ( tfac, dtfac_dt, dt2fac_dt2, &
                                      exact_initial_data_lmax )
    end if
    print*,'Setting particle position'
    call orbit_info%get_orbit ( r, phi, ur, En, Lz )
    call my_eq%effs%set_particle_pos ( r, phi, ur, En, Lz,  &
                               0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp )
    print*,'Particle position set'
    print*,'Evaluating effective source'
    call my_eq%effs%evaluate_source ( my_eq%r_schw, wt )
    print*,'Effective source evaluated'
  end if
  if (use_exact_initial_data) then
    call my_eq%read_all_modes ( )
  end if

  do i = 1, my_eq%effs%nmodes
    call my_eq%effs%source(1,i)%output(rho)
!    if (use_particle) then
!      call my_eq%eq_lcoeffs(i)%output(rho)
!    end if
  end do
  call my_eq%r_schw%output(rho)
!  do i = 1, 4
!    call my_eq%eq_coeffs(i)%output(rho)
!  end do
!  do i = 1, 2
!    call my_eq%eq_lambda(i)%output(rho)
!    do j = 1, 2
!      call my_eq%eq_s(i,j)%output(rho)
!      call my_eq%eq_sinv(i,j)%output(rho)
!    end do
!  end do 

  do i = 1, my_eq%nmodes
    do j = 1, 3
      call my_eq%eq_data(j,i)%output(rho)
      call my_eq%eq_rhs_data(j,i)%output(rho)
    end do
!    do j = 1, 2
!      call my_eq%eq_flux_data(j,i)%output(rho)
!    end do
  end do
  print*,'Calling orbit output'
  call orbit%output()

  allocate ( cobs(3,my_eq%nmodes) )
  do i = 1, my_eq%nmodes
    do j = 1, 3
      call cobs(j,i)%init( (/ rho_min, rho_max /), rho, my_eq%eq_data(j,i) )
      call cobs(j,i)%extract ( )
      call cobs(j,i)%output ( )
    end do
  end do
  print*,'elem_index = ', cobs(1,1)%elem_index
  print*,'node_index = ', cobs(1,1)%node_index

  if (use_particle) then
    call sfobs%init( (/ rho_particle /), rho, my_eq )
    call sfobs%extract ( )
    call sfobs%output ( )
    call sf%output ( )
  end if

  call choose_integrator ( )

  if (use_particle .and. use_generic_orbit) then
      call set_dtime ( coufac*delta_rho_min/my_eq%time_dep_coord%maxspeed )
  else
    call set_dtime ( coufac*delta_rho_min )
  end if
  print*,'dtime = ', get_current_dtime ( )

  call mol_int%init(eqs)
  print*,'ntemp for integrator ', mol_integrator, ' = ', mol_int%ntemp()

  if ( use_particle ) then
    print*,'turn_on_source_smoothly = ', turn_on_source_smoothly
    print*,'exact_initial_data_lmax = ', exact_initial_data_lmax, lmax
    if ( turn_on_source_smoothly .or. &
         (use_exact_initial_data .and. exact_initial_data_lmax < lmax ) ) then
      if (get_current_dtime()*100.0_wp>4.0_wp*tsigma) then
        print*,'In order to turn on the source smoothly, using dtime = ', 4.0_wp*tsigma/100
        call set_dtime ( 4.0_wp*tsigma/100 )
        short_timesteps_active = .true.
      end if
    end if
  end if

  print*,'Starting evolution at time = ', get_current_time ( )
  k = 0
  do 
    if ( get_current_time () > t_final ) exit
    k = k+1
!    print*,'time = ', get_current_time ( )
    call mol_int%step( )

    if (out1d_every>0 .and. mod(k,out1d_every)==0) then
      do i = 1, my_eq%nmodes
        do j = 1, 3
          call my_eq%eq_data(j,i)%output(rho)
          call my_eq%eq_rhs_data(j,i)%output(rho)
        end do
!        do j = 1, 2
!!          call my_eq%eq_flux_data(j,i)%output(rho)
!        end do
        if (use_particle) then
          call my_eq%effs%source(1,i)%output(rho)
        end if
      end do
!      do j = 1, 4
!        call my_eq%eq_coeffs(j)%output(rho)
!      end do 
!      call my_eq%r_star%output(rho)
!      call my_eq%r_schw%output(rho)
    end if
    if (out0d_every>0 .and. mod(k,out0d_every)==0) then
      do i = 1, my_eq%nmodes
        do j = 1, 3
          call cobs(j,i)%extract ( )
          call cobs(j,i)%output ( )
        end do
      end do
      if (use_particle) then
        call sfobs%extract ( )
        call sfobs%output ( )
      end if
      call orbit%output ( )
      call sf%output ( )
    end if
  end do

  print*,'Program ended after ', k, ' iterations.'
end program test