sll_m_linear_pendulum_operators.F90

Go to the documentation of this file.

 
 
module sll_m_linear_pendulum_operators
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#include "sll_memory.h"
#include "sll_working_precision.h"
 
   use sll_m_operator_splitting, only: &
      sll_s_do_split_steps, &
      sll_s_operator_splitting_init, &
      sll_t_operator_splitting
 
   implicit none
 
   public :: &
      sll_s_check_order
 
   private
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
   type, extends(sll_t_operator_splitting) :: linear_pendulum_operators
      sll_real64 :: x  
      sll_real64 :: v  
   contains
      procedure, pass(this) :: operatort => push_x 
      procedure, pass(this) :: operatorv => push_v 
   end type linear_pendulum_operators
 
   ! module global variables
   sll_real64, parameter :: omega = 2.0_f64       
   sll_real64, parameter :: x0 = 1.0_f64          
   sll_real64, parameter :: v0 = 2.0_f64          
   sll_real64, parameter :: t_final = 1.0_f64     
contains
 
   subroutine push_x(this, dt)
      class(linear_pendulum_operators), intent(inout)   :: this
      sll_real64, intent(in)                            :: dt   
 
      this%x = this%x + this%v*dt
   end subroutine push_x
 
   subroutine push_v(this, dt)
      class(linear_pendulum_operators), intent(inout)   :: this
      sll_real64, intent(in)                            :: dt   
 
      this%v = this%v - omega**2*this%x*dt
   end subroutine push_v
 
   subroutine sll_s_check_order(method, steps_fine, expected_order, test_passed)
      sll_int32, intent(in)  :: method         
      sll_real64, intent(in) :: steps_fine     
      sll_int32, intent(in)  :: expected_order 
      logical, intent(inout) :: test_passed
 
      ! local variables
      sll_real64 :: dt
      sll_int32  :: number_time_steps
      sll_real64 :: x_exact, v_exact
      sll_real64 :: error0, error1, error2, order1, order2
      type(linear_pendulum_operators), target :: my_pendulum
      class(sll_t_operator_splitting), pointer :: split
      !sll_int32 :: ierr
 
      ! initialise my_pendulum
      split => my_pendulum
      !  call sll_s_operator_splitting_init(split,sll_p_lie_vt)
      !  call sll_s_operator_splitting_init(split,sll_p_strang_tvt)
      !  call sll_s_operator_splitting_init(split,sll_p_triple_jump_tvt)
      call sll_s_operator_splitting_init(split, method)
 
      ! compute exact solution at final time
      x_exact = x0*cos(omega*t_final) + (v0/omega)*sin(omega*t_final)
      v_exact = -x0*omega*sin(omega*t_final) + v0*cos(omega*t_final)
 
      ! do iterations with smallest time step
      my_pendulum%x = x0
      my_pendulum%v = v0
      dt = t_final/steps_fine
      number_time_steps = int(steps_fine, i32)
 
      call sll_s_do_split_steps(split, dt, number_time_steps)
 
      ! compute  mean square error
      error0 = sqrt((my_pendulum%x - x_exact)**2 + (my_pendulum%v - v_exact)**2)
 
      ! do iterations with middle time step
      my_pendulum%x = x0
      my_pendulum%v = v0
      dt = 2*dt
      number_time_steps = number_time_steps/2
      call sll_s_do_split_steps(split, dt, number_time_steps)
 
      ! compute mean square error
      error1 = sqrt((my_pendulum%x - x_exact)**2 + (my_pendulum%v - v_exact)**2)
 
      ! do iterations with largest time step
      my_pendulum%x = x0
      my_pendulum%v = v0
      dt = 2*dt
      number_time_steps = number_time_steps/2
      call sll_s_do_split_steps(split, dt, number_time_steps)
 
      ! compute  mean square error
      error2 = sqrt((my_pendulum%x - x_exact)**2 + (my_pendulum%v - v_exact)**2)
 
      ! compute order
      !print*, 'error', error0, error1, error2
      order1 = log(error1/error0)/log(2.0_f64)
      order2 = log(error2/error1)/log(2.0_f64)
      if ((abs(order1 - expected_order) > 5.e-2) .or. (abs(order2 - expected_order) > 5.e-2)) then
         test_passed = .false.
         print *, 'error coarse =', error2
         print *, 'error middle =', error1
         print *, '      order (coarse/middle) =', order2
         print *, 'error fine   =', error0
         print *, '      order (middle/fine) =', order1
      end if
   end subroutine sll_s_check_order
 
end module sll_m_linear_pendulum_operators