! Copyright (c) 2024, Paul Muir, Jack Pew, Connor Tannahill
! 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.
! This file contains the BACOLI source code (including the new code
! that implements the interpolation based spatial error estimates.)

module rootfinding
! This module enables the main program (below) to interact with
! the BACOLIKR software in order to implement root finding
! while computing a numerical solution to a PDE. The root finding
! routine that implements the gstop capability employs the 
! bacoli95_vals routine from BACOLIKR.

    use bacoli95_mod, only: bacoli95_sol, bacoli95_vals

    implicit none

!   Global variables and data structures
    integer, parameter          :: dp = kind(0d0)

    ! Structured type containing numerical solution info, accessed 
    ! in this module for convenience 
    type(bacoli95_sol), pointer :: sol

    ! Work array for use in root subroutine, rt, so that calls to
    ! the 'values' subroutine can be employed
    real(dp), pointer           :: work(:)

    contains

    ! 1. 'rt' subroutine. This routine is called in
    ! order to evaluate the gstop vector. This routine is called
    ! by DASKR at the end of every successful time step to 
    ! check if any components of the gstop vector have changed
    ! sign. If a sign change is detected, the 'rt' subroutine is
    ! then called within a search algorithm to locate the time
    ! at which some component of the gstop vector changes sign.

    subroutine rt(t, x, nint, ub, utb, neq, rval, nrt)
    ! t is the current time, x is the spatial mesh, nint is
    ! the number of spatial mesh subintervals, ub is the vector
    ! of B-spline coefficients, ubt is the vector of time
    ! derivatives of the B-spline coefficients, neq is the
    ! number of B-spline coefficieints, rval is the vector
    ! of gstop values at time t, nrt is the number of components
    ! of the gstop vector.

    ! In this example we are interested in determing the time
    ! at which the first solution component reaches steady state. For 
    ! this problem, Gierer-Meinhardt, we will detect a steady state by 
    ! looking for the time at which the first solution has a time 
    ! derivative that is sufficiently close to zero for all x across  
    ! the problem domain.

    ! A practical interpretation of this condition is that the absolute
    ! value of the integral of the time derivative of the first solution 
    ! component over the spatial domain should be equal to 
    ! the user tolerance.

    ! This routine accesses the sol structured type and work array globally

        ! Input arguments
        integer ::  nint, neq, nrt 
        real(dp) :: t, x(nint+1), ub(neq), utb(neq)
       
        ! Output arguments
        real(dp) :: rval(nrt)

        ! Local variables
        real(dp)            :: approx(2)  ! Integral approximation
        real(dp)            :: int        ! Target value of the integral
        real(dp)            :: ut(2,6500) ! Values of time derivative of  
                                          ! solution; assumes nint_max <= 500
        real(dp)            :: h(500)     ! h values; assumes nint_max <= 500
        real(dp)            :: rho(13)    ! Gauss points per subinterval
        real(dp)            :: wts(13)    ! Gauss weights per subinterval
        real(dp)            :: excol(6500), ewts(6500) 
		            ! Total set of Gauss points and weights over 
                            ! entire domain; assumes nint_max <= 500.
                
        integer             :: npts     ! Number of evaluation points
        integer, parameter  :: nderiv = 0 
                               ! Number of spatial derivatives of
                               ! solution that are required
							   
        integer             :: i, ii, j        ! loop indexes

        !-----------------------------------------------------------

	! We are looking for the time when the integral of the absolute
        ! value of the time derivative of the first solution component 
        ! is less than or equal to the user tolerance so the target value 
        ! for the integral is 1d-6. 

        int = 1.0d-6

	! For efficiency, the computation of the Gauss pts and Wts should
	! be done once beforehand and then the array can be available
	! globally. For now we leave these computations here.
        ! For the kcol+3 Gaussian Quadrature rule approximation to the
        ! integral we will need to evaluate the approximate solution
        ! at nint*(kcol+3) points, i.e., the kcol+3 Gauss points on 
	! each subinterval. Therefore

        npts = nint*(sol%kcol+3)
		
	! Next we compute the entire set of Gauss points and weights 
	! for the entire spatial domain. The code below is taken from bacoli.f.
			
	! Calculate the subinterval size sequence.
        do i = 1, nint
            h(i) = x(i+1)-x(i)
        end do

	! Compute the kcol+3 Gaussian points and Gaussian weights on [0,1].
	! (In the call below to the F77 routine, gauleg, the work array
	! needs to be of dimension at least (kcol+3)**2, which has a max
	! of 169 (since the max kcol value is 10). The size of work is
	! set in the next routine below, called setsol. As long as nint_max
	! is large enough, work will be large enough for the call to gauleg.)
	! The points and weights are returned in rho and wts resepctively.
         call gauleg(sol%kcol+3, (sol%kcol+3)*(sol%kcol+3), rho, wts, work, 4)

        ! Define the Gaussian quadrature point and weight sequence over the 
	! whole domain
        do i = 1, nint
            ii = (i - 1) * (sol%kcol + 3)
            do j = 1, sol%kcol+3
                excol(ii + j) = x(i) + h(i) * rho(j)
                ewts(ii + j) = h(i) * wts(j)
            end do
        end do
		
        ! Evaluate the time derivative of the approximate solution to obtain 
        ! values for the integral. It turns out that 'values' can compute all of 
        ! these values in one call.
 
        ! Here we call 'values' with input vector 'utb' which evaluates
        ! the time derivative of the approximate solution at excol and at the 
        ! current time. Replacing 'utb' with 'ub' would evaluate the approximate 
        ! solution at excol and the current time. The 'values' routine is an F77 
        ! routine that can be accessed globally.
        call values(sol%kcol, excol, nint, x, sol%npde, npts, nderiv,  &
                                                             ut, utb, work)

        ! The values of the time derivative of the approximate solution at excol and 
        ! the current time are returned in ut. Inside 'values' the formal argument, 
        ! usol, is declared as follows:
	!        double precision        usol(npde, npts, nderiv+1)
        ! usol contains the solution and the spatial partial derivatives (up to the 
        ! nderiv-th derivative) at x values contained in xsol and at the current time.
        ! Since we have nderiv = 0, usol is effectively an npde by npts array. The
        ! ith column of usol contains the approximate values for the two solution 
        ! components evaluated at the ith xsol value and the current time.

        ! Apply the Gauss rule to obtain the integral approximation
        approx = 0.0d0
        do i = 1, nint
            ii = (i - 1) * (sol%kcol + 3)
            do j = 1, sol%kcol+3
                approx(1:2) = approx(1:2) + ewts(ii+j)*abs(ut(1:2,ii+j))
            end do
        end do
        
        write(7,*) t, approx

        ! We next assign the gstop value to be difference between the
        ! approximate integral value associated with the time derivative of 
        ! the first solution component and the target integral value. 
        rval(1) = approx(1) - int

        ! Initially the time derivative will be positive since the
        ! first solution component initially grows; as the integration 
        ! proceeds the first solution begins to decrease in magnitude 
        ! and spikes begin to form. For the parameter choices we employ
        ! these spikes are stable and the first solution component
        ! converges to a steady state, at which point the time derivative,
        ! across the entire spatial domain becomes close to zero.

    end subroutine rt

    ! 2. Subroutine to allocate work array for use in 'rt' subroutine when it calls
    ! the 'values' routine to evalute the approximate solution. The required 
    ! dimension of the work array depends upon some of the fields of the 'sol' 
    ! structured type so this routine takes 'sol' as an input argument
    subroutine setsol(in_sol)
        type(bacoli95_sol), target :: in_sol

        ! Local variables used in defining size of work array
        integer, parameter :: nconti = 2
        integer, parameter :: nderiv = 0

        ! Assign the pointer in the module scope to be equal to in_sol
        sol => in_sol
        
        ! Allocate work memory for use in 'rt'
        allocate(work((sol%kcol+nconti)*(nderiv+1)+(sol%kcol+nconti+1)   &
            *(sol%kcol+nconti+2)/2+2*sol%kcol*(sol%nint_max+1)+2*nconti))
    end subroutine setsol

end module rootfinding


program simple_example_driver
    ! In this example, all references to BACOLI95 are meant to be
    ! understood as referring to BACOLIKR95. Likewise BACCOLI
    ! references actually refer to BACOLIKR.

    ! A minimal example driver for the BACOLI95 wrapper of BACOLI.
    ! See the comments within bacoli95.f95 for details regarding use.
    use bacoli95_mod, only: bacoli95_init, bacoli95, bacoli95_vals, &
                            bacoli95_sol, bacoli95_sol_teardown

    use rootfinding, only: rt, setsol

    implicit none

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

    type(bacoli95_sol), target  :: sol

    ! Number of PDEs = 2, Number of output points = 101, Max number
    ! of subintervals allowed in the spatial mesh = 500, number of
    ! solution spatial derivatives required, nderiv = 0.
    integer,  parameter :: npde = 2, nout = 1001, nint_max = 500, nderiv = 0

    ! xa is the left spatial boundary, xb is the right spatial boundary
    real(dp), parameter    :: xa = -2.0d0, xb = 2.0d0

    ! Terminal time value = 1000
    real(dp), parameter :: tout = 1000.0d0

    ! xout(:) contains output points, uout(:) contains solution
    ! approximations at output points
    real(dp) :: xout(nout), uout(2,nout)

    ! atol(1) is absolute tolerance, rtol(1) is relative tolerance
    real(dp) :: atol(1), rtol(1)

    ! loop index, i, nrt is the number of gstop roots to be found
    integer :: i, nrt

    ! Declare initial mesh of 20 subintervals
    real(dp)            :: InitMesh(21)

    !   PROBLEM CONSTANTS
    real(dp)            ::   L, D, eps, tau, mu
    COMMON /GMSTUFF/         L, D, eps, tau, mu

!   Subroutines that define the rhs of the PDE, the boundary conditions,
!   the initial condition, and the derivatives of each of these.
    external f, bndxa, bndxb, uinit, derivf, difbxa, difbxb

!----------------------------------------------------------------------

    L = 1.0d0
    D = 0.1d0
    eps = 0.1d0
    mu = 1.0d0
    tau = 1.0d0

!   Number of gstop roots
    nrt = 1

!   Set tolerance
    atol(1) = 1d-6
    rtol(1) = atol(1)

    ! Initialization of spatial mesh
    InitMesh = (/-1.0D0,-0.9D0,-0.8D0,-0.7D0,-0.6D0, -0.5D0,-0.4D0, -0.3D0,-0.2D0, &
                -0.1D0,0.0D0,0.1D0,0.2D0,0.3D0, 0.4D0,0.5D0,0.6D0, 0.7D0,0.8D0,0.9D0, &
                 1.0D0/)
    InitMesh = InitMesh * 2

    ! Initialize the sol structured type
    ! Initial mesh is uniform partition of [-2,2], error estimator is the LOI,
    ! number of collocation points per subinterval, kcol, is 4
    call bacoli95_init(sol, npde, InitMesh, nint_max=nint_max, &
               nrt=nrt, estimator=0, kcol=5, atol=atol, rtol=rtol)

    ! Allocate memory for the work array to be used during the root finding
    call setsol(sol)

    write(6,*) 'THE INPUT IS'
    write(6,*) sol%kcol, sol%nint, sol%npde, tout
    write(6,*) sol%atol(1), sol%rtol(1), "LOI"

    ! Solve the PDE to tout or stop early if a gstop root is found.

    !call bacoli95(sol, tout, f, bndxa, bndxb, uinit, derivf, difbxa, difbxb, rt=rt)
    call bacoli95(sol, tout, f, bndxa, bndxb, uinit, rt=rt)
    ! Stop due to an error condition
    if (sol%idid .le. 0) then
            write(*,*) 'Error? Idid = ', sol%idid
            stop 
    end if

    ! Successful return from bacoli95. We are stopping at a time
    ! that corresponds to finding a gstop root or at tout if no root 
    ! was found.
    if (sol%idid == 5) then
            print'("done. idid=",i5)',sol%idid
            write(*,*) 'Root found at t = ', sol%t0
    else
            print'("done. idid=",i5)',sol%idid
            write(*,*) 'Root not found before reaching tout'
    end if

    ! sol%idid > 0 implies a successful return from bacoli95.
    if (sol%idid > 0) then

        write(*,*) 'Output solution'

        xout = (/xa,(xa+i*(xb-xa)/(nout-1), i=1,nout-1)/)
        call bacoli95_vals(sol, xout, uout)

        write(*,*)

        print'("at t=",f8.2)', sol%t0
        print'(5x,"x",5x5x,"u(t,x)")'
        print '(/"THE OUTPUT IS")'
        print '(a13,a27)', "XOUT", "UOUT"
        do i = 1, nout
            print*, xout(i), uout(:,i)
        end do

    end if

    write(*,*) 'stats'
    write(*,*) 'remeshings: ', sol%num_remeshings
    write(*,*) 'num_ini_remeshings: ', sol%num_ini_remeshings
    write(*,*) 'num_cold_starts: ', sol%num_cold_restarts
    write(*,*) 'num_accepted_time_steps: ', sol%num_accepted_time_steps

    call bacoli95_sol_teardown(sol)

end program

