testSparse.f90 Source File

program testSparse_prg
	use kinds_mod
	use sparse_mod
	use basicSolvers_mod
	use solvers_mod
	use array_mod
	use plplotlib_mod
	implicit none
	
	call setup(fileName='testsSparse-%n.svg',figSize=[400,300])
	call testNewSparse
	call testSpvec
	call testBasicSolvers
	call testSolvers
	call show()
	
contains

	subroutine testNewSparse
		type(sparse_t)::A
		integer::N,M
		
		N = 3
		M = 5
		
		A = newSparse(N,M)
	end subroutine testNewSparse

	subroutine testSpvec
		type(spvec_t)::u,v,r
		
		u%i = [1,2,3]
		u%v = [1.0_wp,2.0_wp,3.0_wp]
		
		v%i = [2,3,4]
		v%v = [2.0_wp,3.0_wp,4.0_wp]
		
		write(*,*) u.o.v
		
		r = u+v
		write(*,*) r%i
		write(*,*) r%v
		
		r = 2.0_wp*u*2.0_wp
		write(*,*) r%i
		write(*,*) r%v
	end subroutine testSpvec

	subroutine testBasicSolvers
		real(wp),parameter::Tl  = 0.0_wp
		real(wp),parameter::Tr  = 1.0_wp
		real(wp),parameter::k   = 1.0_wp
		real(wp),parameter::q0  = 10.0_wp
		real(wp),parameter::tol = 1.0E-8_wp
		
		real(wp)::Ap,Ae,Aw,dx
		type(sparse_t)::A
		real(wp),dimension(:),allocatable::x,q
		real(wp),dimension(:),allocatable::T1,T2,T3
		integer::N,i,s
		
		N  = 1000
		s  = N/30
		allocate(x(0:N+1))
		x  = linspace(0.0_wp,1.0_wp,N+2)
		q  = [( q0     , i=1,N )]
		A  = newSparse(N,N)
		
		allocate(T1(0:N+1))
		allocate(T2(0:N+1))
		allocate(T3(0:N+1))
		
		T1(0)   = Tl
		T1(N+1) = Tr
		
		T2(0)   = Tl
		T2(N+1) = Tr
		
		T3(0)   = Tl
		T3(N+1) = Tr
		
		dx = x(2)-x(1)
		
		i = 1
		Ae = k/dx**2
		Aw = k/dx**2
		Ap = Ae+Aw
		call A%set(i,i  , Ap)
		call A%set(i,i+1,-Ae)
		q(i) = q(i)+Aw*Tl
		do i=2,N-1
			Ae = k/dx**2
			Aw = k/dx**2
			Ap = Ae+Aw
			call A%set(i,i-1,-Aw)
			call A%set(i,i  , Ap)
			call A%set(i,i+1,-Ae)
		end do
		i = N
		Ae = k/dx**2
		Aw = k/dx**2
		Ap = Ae+Aw
		call A%set(i,i-1,-Aw)
		call A%set(i,i  , Ap)
		q(i) = q(i)+Ae*Tr
		
		T1(1:N) = biConjugateGradientStabilized(A,q)
		T2(1:N) = conjugateGradient(A,q)
		T3(1:N) = successiveOverRelaxation(A,q,1.995_wp)
		
		
		call figure()
		call subplot(1,1,1)
		call xylim(mixval(x),mixval(T1)+[0.0_wp,0.05_wp]*span(T1))
		call plot(x(::s),T1(::s),lineStyle='-',lineColor='r',markStyle='x',markColor='r')
		call plot(x(::s),T2(::s),lineStyle='-',lineColor='b',markStyle='o',markColor='b')
		call plot(x(::s),T3(::s),lineStyle='-',lineColor='c',markStyle='s',markColor='c')
		call ticks()
		call labels('Position #fix#fn','Temperature #fiT#fn','1D Heat Conduction with Generation')
		
	end subroutine testBasicSolvers

	subroutine testSolvers
		real(wp),parameter::Tl  = 0.0_wp
		real(wp),parameter::Tr  = 1.0_wp
		real(wp),parameter::k   = 1.0_wp
		real(wp),parameter::q0  = 10.0_wp
		real(wp),parameter::tol = 1.0E-8_wp
		
		real(wp)::Ap,Ae,Aw,dx
		type(sparse_t)::A
		class(solver_t),allocatable::solver
		real(wp),dimension(:),allocatable::x,q
		real(wp),dimension(:),allocatable::T1
		integer::N,i,s
		
		N  = 100
		s  = N/30
		allocate(x(0:N+1))
		x  = linspace(0.0_wp,1.0_wp,N+2)
		q  = [( q0     , i=1,N )]
		A  = newSparse(N,N)
		
		allocate(T1(0:N+1))
		
		T1(0)   = Tl
		T1(N+1) = Tr
		
		dx = x(2)-x(1)
		
		i = 1
		Ae = k/dx**2
		Aw = k/dx**2
		Ap = Ae+Aw
		call A%set(i,i  , Ap)
		call A%set(i,i+1,-Ae)
		q(i) = q(i)+Aw*Tl
		do i=2,N-1
			Ae = k/dx**2
			Aw = k/dx**2
			Ap = Ae+Aw
			call A%set(i,i-1,-Aw)
			call A%set(i,i  , Ap)
			call A%set(i,i+1,-Ae)
		end do
		i = N
		Ae = k/dx**2
		Aw = k/dx**2
		Ap = Ae+Aw
		call A%set(i,i-1,-Aw)
		call A%set(i,i  , Ap)
		q(i) = q(i)+Ae*Tr
		
		allocate(solver,source=jacobi_t())
		call solver%setup(A)
		T1(1:N) = solver%solve(A,q)
		deallocate(solver)
		
		allocate(solver,source=gaussSeidel_t())
		call solver%setup(A)
		T1(1:N) = solver%solve(A,q)
		deallocate(solver)
		
		allocate(solver,source=SOR_t(1.99_wp))
		call solver%setup(A)
		T1(1:N) = solver%solve(A,q)
		deallocate(solver)
		
		allocate(solver,source=conjugateGradient_t())
		call solver%setup(A)
		T1(1:N) = solver%solve(A,q)
		deallocate(solver)
		
		allocate(solver,source=biCGSTAB_t())
		call solver%setup(A)
		T1(1:N) = solver%solve(A,q)
		deallocate(solver)
		
		call figure()
		call subplot(1,1,1)
		call xylim(mixval(x),mixval(T1)+[0.0_wp,0.05_wp]*span(T1))
		call plot(x(::s),T1(::s),lineStyle='-',lineColor='k',markStyle='x',markColor='r')
		call ticks()
		call labels('Position #fix#fn','Temperature #fiT#fn','1D Heat Conduction with Generation')
		
	end subroutine testSolvers

end program testSparse_prg