module dpres
  integer, parameter :: dp=kind(1.D0)
end module dpres
!Subroutine calculant l equation a resoudre par differences finies
!A.u = [Id - k d2/dx2].u + CL
subroutine equation(N,dx,k,u,Au)
use dpres
implicit none
integer :: N,i
real(kind=dp) :: d2udx2,dx
real(kind=dp), dimension(0:N-1) :: k,u,Au

do i=1,N-2
  d2udx2 = (u(i+1) - 2.D0*u(i) + u(i-1))/(dx*dx) ! calcul de derivees seconde par DF
  Au(i) = u(i)-k(i)*d2udx2 ! calcul de l equation
enddo
!Application des CL
Au(0) = u(0) ! condition de dirichlet homogene
Au(N-1) = (u(N-1) - u(N-2))/dx; ! condition de Neumann homogène
end subroutine equation
!subroutine de calcul d un produit scalaire
subroutine pscal(N,u,v,somme)
use dpres
implicit none
integer :: i,N
real(kind=dp) :: somme
real(kind=dp), dimension(0:N-1) :: u,v
somme = 0.D0
do i=0,N-1
  somme = somme + u(i)*v(i)
enddo
end subroutine pscal
!subroutine de calcul du max d'un tableau
subroutine maxi(N,u,maximum)
use dpres
implicit none
integer :: i,N
real(kind = dp) :: maximum
real(kind=dp), dimension(0:N-1) :: u
maximum = dabs(u(0))
do i=1,N-1
  if(dabs(u(i))>maximum) then
    maximum = dabs(u(i))
  endif
enddo
end subroutine maxi
!schema iteratif de Richardson
!par minimisation du residu
subroutine richardson(N,dx,k,uj,smb,itmax)
use dpres
implicit none
integer :: i,j=0,N,itmax
real(kind=dp) :: dx,alpha,numer,denom,res
real(kind=dp), dimension(0:N-1) :: k,uj,smb,rj,Arj ! rj = erreur point a point, Arj = L(rj)
!calcul du residu initial
call equation(N,dx,k,uj,rj)
do i=0,N-1
  rj(i) = smb(i) - rj(i) ! rj = smb - A.uj
enddo
call maxi(N,rj,res)
do while(res.gt.1.d-6.and.j<itmax)
  call equation(N,dx,k,rj,Arj) ! On applique l equation sur l erreur
  call pscal(N,rj,Arj,numer) !<rj,A.rj>
  call pscal(N,Arj,Arj,denom)!<A.rj,A.rj>
  alpha = numer/denom ! facteur de normalisation <rj,A.rj>/<A.rj,A.rj>
  do i = 0,N-1 ! calcul de l itere suivants
    uj(i) = uj(i) + alpha*rj(i) ! iteration sur u
    rj(i) = rj(i) - alpha*Arj(i)! mise a jour de erreur
  enddo
  call maxi(N,rj,res) ! nouveau residu
  j = j + 1
enddo
write(*,*) "j = ",j," r = ",res
end subroutine richardson
!subroutine d ecriture dans un fichier : x(i) et u(i)
subroutine ecriture(ifile,N,x,u)
use dpres
implicit none
integer :: i,N,ifile
real(kind=dp), dimension(0:N-1) :: x,u
do i=0,N-1
  write(ifile,*) x(i),u(i)
enddo
end subroutine ecriture

!programme principal
program thermNL
use dpres
implicit none
integer :: N,Nt,i
real(kind=dp) :: dx,pi=4.d0*datan(1.D0),L = 2.d0*datan(1.D0)
real(kind=dp), allocatable, dimension(:) :: u,k,f,x
!lecture des donnes
open(unit=11,file="param.dat",status="old")
read(11,*) N
read(11,*) Nt
write(*,*) "N = ",N," Nt = ",Nt
!allocation dynamique
allocate(u(0:N-1),x(0:N-1),f(0:N-1),k(0:N-1))
! initialisation des parametres et des champs
dx = L/dfloat(N-1)
do i=0,N-1
  x(i) = dfloat(i)*dx ! espace
  f(i) = 0.d0         ! second membre vide
  k(i) = (1.d0+2.d-1*dexp(-x(i)*x(i)/(0.5d0*0.5d0))) ! terme de raideur NL
  u(i) = cos(2.D0*x(i))  ! initialisation pour calcul du second membre
enddo
!calcul du second membre
call equation(N,dx,k,u,f)
f(0) = 1.d0 ! CL homogene aux deux bords
f(N-1) = 0.d0

do i=0,N-1 ! reinitialisation de la condition initiale
  u(i) = 1.d0
enddo

call richardson(N,dx,k,u,f,Nt) ! resolution iterative
open(unit=12, file="solf.dat",status="unknown")
call ecriture(12,N,x,u) ! sauvegarde
close(11) ; close(12)
deallocate(x,u,f,k)

end program thermNL
