! Driver for EEE2022-Ex3

! Copyright (c) 2013, Paul Muir, Jack Pew
! Paul Muir, Mathematics and Computing Science, Saint Mary's University.
! All rights reserved.
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions
! are met:
! * Redistributions of source code must retain the above copyright
!   notice, this list of conditions and the following disclaimer.
! * Redistributions in binary form must reproduce the above copyright
!   notice, this list of conditions and the following disclaimer in the
!   documentation and/or other materials provided with the distribution.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
! "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
! A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
! HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
! SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
! LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
! DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
! THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
! (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
! OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

program simple_example_driver

    ! A minimal example driver for the EBACOLI95 wrapper of EBACOLI.
    ! See the comments within ebacoli95.f95 for details regarding use.

    use ebacoli95_mod, only: ebacoli95_init, ebacoli95, ebacoli95_vals, &
                            ebacoli95_sol, ebacoli95_sol_teardown

    implicit none
    integer, parameter  :: dp = kind(0d0)

    ! Declare "sol" data structure that will hold solution
    type(ebacoli95_sol)  :: sol

    ! Set "nout", number of output points over spatial domain
    ! and "nint_max", maximum number of subintervals.
    integer,  parameter :: nout = 11, nint_max = 500

    ! Set "npde_sub" to indicate number of pdes, time-dependent ODEs, and
    ! space-dependent ODEs, respectively.
    integer,  parameter, dimension(3) :: npde_sub = (/1,0,1/) !<=

    ! Set "npde" to indicate total number of differential equations 
    integer,  parameter :: npde = 2 !<=

    ! Set "xa" and "xb", the endpoints of spatial domain
    real(dp), parameter :: xa = 0, xb = 1 !<=

    ! Initial mesh of "nin" points
    integer,  parameter :: nin = 11
    real(dp) :: xin(nin)

    ! Set output time, "tout"
    real(dp), parameter :: tout = 1.0D0 !<=

    ! Set tolerance, "tol"
    real(dp), parameter :: tol = 1d-6

    ! Declare arrays for output points, "xout", output solution values, "uout",
    ! exact solution values, "exact", and the error, "error".
    real(dp) :: xout(nout), uout(npde*nout), exact(nout), error(nout), err_tot_sq

    ! Loop variables
    integer :: i, k, ij

    external f, bndxa, bndxb, uinit, truu

    ! Initialize a uniform grid to pass to ebacoli
    xin = (/(xa+i*(xb-xa)/(nin-1), i=0,nin-1)/)

    ! Initialize computation
    call ebacoli95_init(sol, npde_sub, xin, atol=(/tol/), &
         rtol=(/tol/), nint_max = nint_max)

    ! Write out sizes of problem components
    call write_subsystem_sizes(sol%npde_sub)

    ! Compute solution at tout
    call ebacoli95(sol, tout, f, bndxa, bndxb, uinit)

    ! Output idid to check for a successful computation
    print '("idid=",i5)',sol%idid
    if (sol%idid > 0) print '("idid > 0 => Successful computation")'

    ! Output solution at tout for nout values of x uniformly
    ! distributed over spatial domain
    if (sol%idid > 0) then
        xout = (/(xa+i*(xb-xa)/(nout-1), i=0,nout-1)/)

        call ebacoli95_vals(sol, xout, uout)

        print '("At t = ",f7.2)', sol%t0
        write(*,'(/a)') 'the solution is'
        write(*,'(a13,a27)') 'XOUT', 'UOUT'
        do i = 1, nout
           ij = (i-1)*npde
           write(*,*) xout(i), (uout(ij+k), k = 1, npde)
        end do
        
        ! Compute and write out discrete L2-norm error over
        ! all solution components at tout
        call write_discrete_L2_error(sol,truu,nout)

        ! Compute discrete L2 norm for u(x,t)
        write(*,*)
        print '("tol = ",e7.2)', tol
        err_tot_sq = 0.0d0

        do i = 1, nout
           call truu(sol%t0, xout(i), exact(i), npde)
           ij = (i-1)*npde
           error(i) = exact(i) - uout(ij+1)
           err_tot_sq = err_tot_sq + error(i)**2
        end do
        print '("Discrete L2 norm of u")'
        write(*,*) sqrt(err_tot_sq)/sqrt(dble(nout))

     end if

    call ebacoli95_sol_teardown(sol)

end program