integrate Subroutine

private subroutine integrate(scheme, a, b, U0, final_time, Dt, iterations, stages, is_fast, solution, error, last_step)

Integrate domain by means of the given scheme.

Arguments

Type	Intent	Attributes		Name
character(len=*),	intent(in)		::	scheme	Selected scheme.
real(kind=R_P),	intent(in)		::	a	a coefficient.
real(kind=R_P),	intent(in)		::	b	b coefficient.
real(kind=R_P),	intent(in)		::	U0	Initial state.
real(kind=R_P),	intent(in)		::	final_time	Final integration time.
real(kind=R_P),	intent(in)		::	Dt	Time step.
integer(kind=I_P),	intent(in)		::	iterations	Number of fixed point iterations.
integer(kind=I_P),	intent(in)		::	stages	Number of stages.
logical,	intent(in)		::	is_fast	Activate fast mode integration.
real(kind=R_P),	intent(out),	allocatable	::	solution(:,:)	Solution at each time step, X-Y.
real(kind=R_P),	intent(out)		::	error	Error (norm L2) with respect the exact solution.
integer(kind=I_P),	intent(out)		::	last_step	Last time step computed.

Source Code

   subroutine integrate(scheme, a, b, U0, final_time, Dt, iterations, stages, is_fast, solution, error, last_step)
   !< Integrate domain by means of the given scheme.
   character(*),           intent(in)    :: scheme        !< Selected scheme.
   real(R_P),              intent(in)    :: a             !< *a* coefficient.
   real(R_P),              intent(in)    :: b             !< *b* coefficient.
   real(R_P),              intent(in)    :: U0            !< Initial state.
   real(R_P),              intent(in)    :: final_time    !< Final integration time.
   real(R_P),              intent(in)    :: Dt            !< Time step.
   integer(I_P),           intent(in)    :: iterations    !< Number of fixed point iterations.
   integer(I_P),           intent(in)    :: stages        !< Number of stages.
   logical,                intent(in)    :: is_fast       !< Activate fast mode integration.
   real(R_P), allocatable, intent(out)   :: solution(:,:) !< Solution at each time step, X-Y.
   real(R_P),              intent(out)   :: error         !< Error (norm L2) with respect the exact solution.
   integer(I_P),           intent(out)   :: last_step     !< Last time step computed.
   class(integrator_object), allocatable :: integrator    !< The integrator.
   type(integrand_lcce)                  :: domain        !< Linear constant coefficients equation field.
   type(integrand_lcce), allocatable     :: previous(:)   !< Previous time steps solutions.
   type(integrand_lcce), allocatable     :: stage(:)      !< Runge-Kutta stages.
   type(integrand_lcce)                  :: buffer        !< Buffer field.
   type(integrand_lcce)                  :: filter        !< Filter displacement.
   integer                               :: step          !< Time steps counter.
   integer                               :: step_offset   !< Time steps counter offset for slicing previous data array.
   real(R_P), allocatable                :: Dts(:)        !< Time steps for variable stepsize.
   real(R_P)                             :: Dt_a          !< Adaptive time step.

   call domain%initialize(a=a, b=b, U0=U0)

   if (allocated(solution)) deallocate(solution) ; allocate(solution(0:domain%integrand_dimension(), 0:int(final_time/Dt)))
   solution = 0.0_R_P
   solution(1:, 0) = domain%output()

   call foodie_integrator_factory(scheme=scheme, integrator=integrator, stages=stages, tolerance=1e2_R_P)
   if (is_fast) call check_scheme_has_fast_mode(scheme=trim(adjustl(scheme)), integrator=integrator)

   if (integrator%is_multistep()) allocate(previous(1:integrator%steps_number()))
   if (integrator%is_multistage()) allocate(stage(1:integrator%stages_number()))

   select type(integrator)
   type is(integrator_adams_bashforth)
      step = 0
      do while(solution(0, step) < final_time .and. step < ubound(solution, dim=2))
        step = step + 1
        if (integrator%steps_number() >= step) then
           domain = [domain%exact_solution(t=step * Dt)]
           previous(step) = domain
        else
           if (is_fast) then
              call integrator%integrate_fast(U=domain,          &
                                             previous=previous, &
                                             buffer=buffer,     &
                                             Dt=Dt,             &
                                             t=solution(0, step-integrator%steps_number():step-1))
           else
              call integrator%integrate(U=domain,          &
                                        previous=previous, &
                                        Dt=Dt,             &
                                        t=solution(0, step-integrator%steps_number():step-1))
           endif
        endif
        solution(0, step) = step * Dt
        solution(1:, step) = domain%output()
      enddo
      last_step = step

   type is(integrator_adams_bashforth_moulton)
      step = 0
      do while(solution(0, step) < final_time .and. step < ubound(solution, dim=2))
        step = step + 1
        if (integrator%steps_number() >= step) then
           domain = [domain%exact_solution(t=step * Dt)]
           previous(step) = domain
        else
          if (is_fast) then
             call integrator%integrate_fast(U=domain,          &
                                            previous=previous, &
                                            buffer=buffer,     &
                                            Dt=Dt,             &
                                            t=solution(0, step-integrator%steps_number():step-1))
          else
             call integrator%integrate(U=domain,          &
                                       previous=previous, &
                                       Dt=Dt,             &
                                       t=solution(0, step-integrator%steps_number():step-1))
          endif
        endif
        solution(0, step) = step * Dt
        solution(1:, step) = domain%output()
      enddo
      last_step = step

   type is(integrator_adams_moulton)
      if (allocated(previous)) deallocate(previous) ; allocate(previous(1:integrator%steps_number()+1))
      if (integrator%steps_number()==0) then
        step_offset = 1                         ! for 0 step-(a convention)-solver offset is 1
      else
        step_offset = integrator%steps_number() ! for >0 step-solver offset is steps
      endif
      step = 0
      do while(solution(0, step) < final_time .and. step < ubound(solution, dim=2))
        step = step + 1
        if (integrator%steps_number() >= step) then
           domain = [domain%exact_solution(t=step * Dt)]
           previous(step) = domain
        else
          if (is_fast) then
             if (iterations>1) then
               call integrator%integrate_fast(U=domain,                              &
                                              previous=previous,                     &
                                              buffer=buffer,                         &
                                              Dt=Dt,                                 &
                                              t=solution(0,step-step_offset:step-1), &
                                              iterations=iterations)
             else
               call integrator%integrate_fast(U=domain,          &
                                              previous=previous, &
                                              buffer=buffer,     &
                                              Dt=Dt,             &
                                              t=solution(0,step-step_offset:step-1))
             endif
          else
             if (iterations>1) then
               call integrator%integrate(U=domain,                              &
                                         previous=previous,                     &
                                         Dt=Dt,                                 &
                                         t=solution(0,step-step_offset:step-1), &
                                         iterations=iterations)
             else
               call integrator%integrate(U=domain,          &
                                         previous=previous, &
                                         Dt=Dt,             &
                                         t=solution(0,step-step_offset:step-1))
             endif
          endif
        endif
        solution(0, step) = step * Dt
        solution(1:, step) = domain%output()
      enddo
      last_step = step

   type is(integrator_back_df)
      step = 0
      do while(solution(0, step) < final_time .and. step < ubound(solution, dim=2))
        step = step + 1
        if (integrator%steps_number() >= step) then
           domain = [domain%exact_solution(t=step * Dt)]
           previous(step) = domain
        else
          if (is_fast) then
             call integrator%integrate_fast(U=domain,          &
                                            previous=previous, &
                                            buffer=buffer,     &
                                            Dt=Dt,             &
                                            t=solution(0, step-integrator%steps_number():step-1))
          else
             call integrator%integrate(U=domain,          &
                                       previous=previous, &
                                       Dt=Dt,             &
                                       t=solution(0, step-integrator%steps_number():step-1))
          endif
        endif
        solution(0, step) = step * Dt
        solution(1:, step) = domain%output()
      enddo
      last_step = step

   type is(integrator_euler_explicit)
      step = 0
      do while(solution(0, step) < final_time .and. step < ubound(solution, dim=2))
        step = step + 1
        if (is_fast) then
           call integrator%integrate_fast(U=domain, buffer=buffer, Dt=Dt, t=solution(0, step))
        else
           call integrator%integrate(U=domain, Dt=Dt, t=solution(0, step))
        endif
        solution(0, step) = step * Dt
        solution(1:, step) = domain%output()
      enddo
      last_step = step
   type is(integrator_leapfrog)
      step = 0
      do while(solution(0, step) < final_time .and. step < ubound(solution, dim=2))
        step = step + 1
        if (integrator%steps_number() >= step) then
           domain = [domain%exact_solution(t=step * Dt)]
           previous(step) = domain
        else
           if (index(scheme, 'raw') > 0 ) then
              if (is_fast) then
                call integrator%integrate_fast(U=domain, previous=previous, buffer=buffer, Dt=Dt, t=solution(0, step),filter=filter)
              else
                call integrator%integrate(U=domain, previous=previous, Dt=Dt, t=solution(0, step), filter=filter)
              endif
           else
              if (is_fast) then
                call integrator%integrate_fast(U=domain, previous=previous, buffer=buffer, Dt=Dt, t=solution(0, step))
              else
                call integrator%integrate(U=domain, previous=previous, Dt=Dt, t=solution(0, step))
              endif
           endif
        endif
        solution(0, step) = step * Dt
        solution(1:, step) = domain%output()
      enddo
      last_step = step

   type is(integrator_lmm_ssp)
      step = 0
      do while(solution(0, step) < final_time .and. step < ubound(solution, dim=2))
        step = step + 1
        if (integrator%steps_number() >= step) then
           domain = [domain%exact_solution(t=step * Dt)]
           previous(step) = domain
        else
          if (is_fast) then
             call integrator%integrate_fast(U=domain,          &
                                            previous=previous, &
                                            buffer=buffer,     &
                                            Dt=Dt,             &
                                            t=solution(0, step-integrator%steps_number():step-1))
          else
             call integrator%integrate(U=domain,          &
                                       previous=previous, &
                                       Dt=Dt,             &
                                       t=solution(0, step-integrator%steps_number():step-1))
          endif
        endif
        solution(0, step) = step * Dt
        solution(1:, step) = domain%output()
      enddo
      last_step = step

   type is(integrator_lmm_ssp_vss)
      if (allocated(Dts)) deallocate(Dts) ; allocate(Dts(1:integrator%steps_number())) ; Dts = Dt
      step = 0
      do while(solution(0, step) < final_time .and. step < ubound(solution, dim=2))
        step = step + 1
        if (integrator%steps_number() >= step) then
           domain = [domain%exact_solution(t=step * Dt)]
           previous(step) = domain
        else
          if (is_fast) then
             call integrator%integrate_fast(U=domain,          &
                                            previous=previous, &
                                            buffer=buffer,     &
                                            Dt=Dts,            &
                                            t=solution(0, step-integrator%steps_number():step-1))
          else
             call integrator%integrate(U=domain,          &
                                       previous=previous, &
                                       Dt=Dts,            &
                                       t=solution(0, step-integrator%steps_number():step-1))
          endif
        endif
        solution(0, step) = step * Dt
        solution(1:, step) = domain%output()
      enddo
      last_step = step

   type is(integrator_ms_runge_kutta_ssp)
      step = 0
      do while(solution(0, step) < final_time .and. step < ubound(solution, dim=2))
        step = step + 1
        if (integrator%steps_number() >= step) then
           domain = [domain%exact_solution(t=step * Dt)]
           previous(step) = domain
        else
          if (is_fast) then
             call integrator%integrate_fast(U=domain,          &
                                            previous=previous, &
                                            stage=stage,       &
                                            buffer=buffer,     &
                                            Dt=Dt,             &
                                            t=solution(0, step-integrator%steps_number():step-1))
          else
             call integrator%integrate(U=domain,          &
                                       previous=previous, &
                                       stage=stage,       &
                                       Dt=Dt,             &
                                       t=solution(0, step-integrator%steps_number():step-1))
          endif
        endif
        solution(0, step) = step * Dt
        solution(1:, step) = domain%output()
      enddo
      last_step = step

   type is(integrator_runge_kutta_emd)
       step = 0
      do while(solution(0, step) < final_time .and. step < ubound(solution, dim=2))
        step = step + 1
        Dt_a = Dt
        if (is_fast) then
           call integrator%integrate_fast(U=domain, stage=stage, buffer=buffer, Dt=Dt_a, t=solution(0, step))
        else
           call integrator%integrate(U=domain, stage=stage, Dt=Dt_a, t=solution(0, step))
        endif
        solution(0, step) = step * Dt
        solution(1:, step) = domain%output()
      enddo
      last_step = step

   type is(integrator_runge_kutta_ls)
      if (allocated(stage)) deallocate(stage) ; allocate(stage(1:integrator%registers_number()))
      step = 0
      do while(solution(0, step) < final_time .and. step < ubound(solution, dim=2))
        step = step + 1
        if (is_fast) then
           call integrator%integrate_fast(U=domain, stage=stage, buffer=buffer, Dt=Dt, t=solution(0, step))
        else
           call integrator%integrate(U=domain, stage=stage, Dt=Dt, t=solution(0, step))
        endif
        solution(0, step) = step * Dt
        solution(1:, step) = domain%output()
      enddo
      last_step = step

   type is(integrator_runge_kutta_lssp)
      step = 0
      do while(solution(0, step) < final_time .and. step < ubound(solution, dim=2))
        step = step + 1
        if (is_fast) then
           call integrator%integrate_fast(U=domain, stage=stage, buffer=buffer, Dt=Dt, t=solution(0, step))
        else
           call integrator%integrate(U=domain, stage=stage, Dt=Dt, t=solution(0, step))
        endif
        solution(0, step) = step * Dt
        solution(1:, step) = domain%output()
      enddo
      last_step = step

   type is(integrator_runge_kutta_ssp)
      step = 0
      do while(solution(0, step) < final_time .and. step < ubound(solution, dim=2))
        step = step + 1
        if (is_fast) then
           call integrator%integrate_fast(U=domain, stage=stage, buffer=buffer, Dt=Dt, t=solution(0, step))
        else
           call integrator%integrate(U=domain, stage=stage, Dt=Dt, t=solution(0, step))
        endif
        solution(0, step) = step * Dt
        solution(1:, step) = domain%output()
      enddo
      last_step = step
   endselect

   error = 0._R_P
   do step=0, last_step
      domain = solution(1:, step)
      error = error + (domain%output() - domain%exact_solution(t=solution(0, step))) ** 2
   enddo
   error = sqrt(error)
   endsubroutine integrate

All Procedures

a a a a a add_burgers add_euler add_euler add_euler add_lorenz allocate_integrand_members allocate_integrand_members allocate_integrand_members allocate_integrand_members allocate_integrand_members amplitude_phase assign_abstract assign_integrand assign_integrand assign_integrand assign_multistage assign_multistage_multistep assign_multistep assign_real assign_real assign_real average_solution burgers_assign_burgers burgers_assign_real burgers_local_error burgers_multiply_burgers burgers_multiply_real check_error check_scheme_has_fast_mode check_scheme_has_fast_mode check_scheme_has_fast_mode class_name class_name class_name class_name class_name class_name class_name class_name class_name class_name class_name class_name class_name compute_dt compute_dt compute_dt compute_dt compute_dt compute_dt compute_dt compute_dx compute_inter_states compute_inter_states compute_inter_states compute_inter_states compute_inter_states conservative2primitive conservative2primitive conservative2primitive conservative2primitive conservative2primitive d2Burgers_dx2 dBurgers_dt dBurgers_dx description description description description destroy destroy destroy destroy destroy destroy destroy destroy destroy destroy destroy destroy destroy destroy destroy destroy destroy destroy destroy destroy_abstract destroy_multistage destroy_multistage_multistep destroy_multistep destroy_rk destroy_rk dEuler_dt dEuler_dt dEuler_dt dEuler_dt dEuler_dt dLorenz_dt dt_ratio dU_dt dU_dt dU_dt E E E E E error error error euler_assign_euler euler_assign_euler euler_assign_euler euler_assign_euler euler_assign_euler euler_assign_real euler_assign_real euler_assign_real euler_local_error euler_local_error euler_local_error euler_multiply_euler euler_multiply_euler euler_multiply_euler euler_multiply_real euler_multiply_real euler_multiply_real exact_solution exact_solution exact_solution execute execute execute export_tecplot export_tecplot export_tecplot foodie_integrator_class_names foodie_integrator_factory foodie_integrator_schemes H H H H H has_fast_mode has_fast_mode has_fast_mode has_fast_mode has_fast_mode has_fast_mode has_fast_mode has_fast_mode has_fast_mode has_fast_mode has_fast_mode has_fast_mode has_fast_mode impose_boundary_conditions impose_boundary_conditions impose_boundary_conditions impose_boundary_conditions impose_boundary_conditions impose_boundary_conditions init init init init init init init init init init init init init init_rk init_rk initialize initialize initialize initialize initialize initialize initialize initialize initialize initialize initialize initialize initialize initialize initialize initialize initialize initialize initialize initialize_order_s initialize_order_s_1 integr_assign_integr integr_assign_integr integr_assign_integr integr_assign_integr integr_assign_integr integr_assign_integr integr_assign_integr integr_assign_integr integr_assign_integr integr_assign_integr integr_assign_integr integr_assign_integr integr_assign_integr integrand_add_integrand integrand_add_integrand integrand_add_integrand integrand_add_integrand_fast integrand_add_integrand_fast integrand_add_integrand_fast integrand_add_real integrand_add_real integrand_add_real integrand_dimension integrand_dimension integrand_dimension integrand_multiply_integrand integrand_multiply_integrand integrand_multiply_integrand integrand_multiply_integrand_fast integrand_multiply_integrand_fast integrand_multiply_integrand_fast integrand_multiply_real integrand_multiply_real integrand_multiply_real integrand_multiply_real_scalar integrand_multiply_real_scalar integrand_multiply_real_scalar integrand_multiply_real_scalar_fast integrand_multiply_real_scalar_fast integrand_multiply_real_scalar_fast integrand_sub_integrand integrand_sub_integrand integrand_sub_integrand integrand_sub_real integrand_sub_real integrand_sub_real integrand_subtract_integrand_fast integrand_subtract_integrand_fast integrand_subtract_integrand_fast integrate integrate integrate integrate integrate integrate integrate integrate integrate integrate integrate integrate integrate integrate integrate integrate integrate_fast integrate_fast integrate_fast integrate_fast integrate_fast integrate_fast integrate_fast integrate_fast integrate_fast integrate_fast integrate_fast integrate_fast integrate_fast integrate_order_2 integrate_order_2_fast integrate_order_3 integrate_order_3_fast integrate_order_s integrate_order_s_1 integrate_order_s_1_fast integrate_order_s_fast integrate_rk integrate_rk is_admissible is_available is_class_available is_multistage is_multistage is_multistage is_multistep is_multistep is_multistep is_scheme_available is_supported is_supported is_supported is_supported is_supported is_supported is_supported is_supported is_supported is_supported is_supported is_supported is_supported local_error local_error local_error lorenz_assign_lorenz lorenz_assign_real lorenz_local_error lorenz_multiply_lorenz lorenz_multiply_real new_Dt observed_order observed_order observed_order omega output output output output output output output output output output p p p p p parse_cli parse_cli parse_cli previous_step previous_step previous_step primitive2conservative primitive2conservative primitive2conservative primitive2conservative primitive2conservative r r r r r real_add_integrand real_add_integrand real_add_integrand real_multiply_burgers real_multiply_euler real_multiply_euler real_multiply_euler real_multiply_integrand real_multiply_integrand real_multiply_integrand real_multiply_lorenz real_scalar_multiply_integrand real_scalar_multiply_integrand real_scalar_multiply_integrand real_sub_integrand real_sub_integrand real_sub_integrand reconstruct_interfaces reconstruct_interfaces_states reconstruct_interfaces_states reconstruct_interfaces_states reconstruct_interfaces_states reconstruct_interfaces_states registers_number riemann_solver riemann_solver riemann_solver riemann_solver riemann_solver save_results save_results save_results save_results save_results save_results save_results save_results save_results save_time_serie save_time_serie save_time_serie save_time_serie save_time_serie scheme_number set_cli set_cli set_cli set_sin_wave_initial_state set_square_wave_initial_state stages_number stages_number stages_number steps_number steps_number steps_number sub_burgers sub_euler sub_euler sub_euler sub_lorenz supported_schemes supported_schemes supported_schemes supported_schemes supported_schemes supported_schemes supported_schemes supported_schemes supported_schemes supported_schemes supported_schemes supported_schemes supported_schemes synchronize synchronize synchronize synchronize t_fast t_fast t_fast test test test_ab test_ab test_ab test_euler test_euler test_euler test_leapfrog test_leapfrog test_leapfrog test_ls_rk test_ls_rk test_ls_rk test_tvd_rk test_tvd_rk test_tvd_rk tokenize trigger_error update_previous update_previous update_previous_steps update_previous_steps update_previous_steps

integrate Subroutine

Source Code

All Procedures

private subroutine integrate(scheme, a, b, U0, final_time, Dt, iterations, stages, is_fast, solution, error, last_step)

Arguments

Calls

Called By

Source Code

Source Code

All Procedures