CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Wiki > Sample code for BiCGSTAB - Fortran 90

Sample code for BiCGSTAB - Fortran 90

From CFD-Wiki

Revision as of 11:02, 27 May 2015 by Alisha (Talk | contribs)
Jump to: navigation, search
	!-----------------------------------MIT LICENCE-----------------------------------------!
	!											!
	!	Copyright (c) 2015, Biswajit Ghosh						!
	!	http://www.biswa.me/cfd								!
	!	All rights reserved.								!
	!											!
	!	Permission is hereby granted, free of charge, to any person obtaining a copy	!
	!	of this software and associated documentation files (the "Software"), to deal	!
	!	in the Software without restriction, including without limitation the rights	!
	!	to use, copy, modify, merge, publish, distribute, sublicense, and/or sell	!
	!	copies of the Software, and to permit persons to whom the Software is		!
	!	furnished to do so, subject to the following conditions:                        !
	!											!
	!	The above copyright notice and this permission notice shall be included in	!
	!	all copies or substantial portions of the Software.				!
	!											!
	!	THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR	!
	!	IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,	!
	!	FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE	!
	!	AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER		!
	!	LIABILITY, WHE  THER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,	!
	!	OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN	!
	!	THE SOFTWARE.									!
	!---------------------------------------------------------------------------------------!

	!---------------------------------------------------------------------------------------!
	!	NOTE : THIS CODE IS DESIGNED ONLY FOR LEARNING PURPOSE.                         !
	!	IMPLEMENTATION IN PRACTICAL PROBLEM MAY RESULT CONVERGENCE ISSUE		!
	!---------------------------------------------------------------------------------------!
	

module var
implicit none
	
	integer, parameter		  	:: n = 4		!-------> size of the matrices
	double precision, parameter		:: e = 1e-20	!-------> maximum permiable error
	double precision, dimension(1:n,1:n)	:: A
	double precision, dimension(1:n    )	:: b,x

	double precision, dimension(1:n    ) 	:: r, rs, v, p, s, t
	double precision		  	:: rho, rho_prev, alpha, omega, beta
	double precision			:: norm_r, norm_b	
end module var

	!---------------------------------------------------------------------------------------!
	!	this code is based on the algo given in wikipedia :				!
	!	http://en.wikipedia.org/wiki/Biconjugate_gradient_stabilized_method		!
	!---------------------------------------------------------------------------------------!


!--------------------------------------main program block---------------------------------------!
program BiCGSTAB
use var
implicit none

	integer 		 :: i,j,it
	double precision :: summesion, temp

	call ab_def	!---------> defination of matrics A,b
	
	call init	!---------> Setting initial values of variable

	it = 0		!---------> with 'it', we count the number of iteration



!*******************************************************!	
	do while(norm_r .GT. e*norm_b)			!-----> convergence criteria
!-------------------------------------------------------! step 01

	rho_prev = rho			!---------> rho_prev stores the 
	rho   = 0.0d0			!			value of rho at n-1 th step
	do i = 1,n
	rho = rho + rs(i)*r(i)
	end do
	
!------------------------------------------------------! step 02

	beta = (rho/rho_prev) * (alpha/omega)
	
!------------------------------------------------------! step 03

	p = r + beta * (p - omega*v)
	
!------------------------------------------------------! step 04

	do i = 1,n
	summesion = 0.0d0
	do j = 1,n
	summesion = summesion + A(i,j)*p(j)
	end do
	v(i)  = summesion
	end do
	
!------------------------------------------------------! step 05

	summesion = 0.0d0
	do i = 1,n
	summesion = summesion + rs(i)*v(i)
	end do
	
	alpha = rho/summesion
	
!------------------------------------------------------! step 06

	s = r - alpha*v
	
!------------------------------------------------------! step 07

	do i = 1,n
	summesion = 0.0d0
	do j = 1,n
	summesion = summesion + A(i,j)*s(j)
	end do
	t(i)  = summesion
	end do
	
!-----------------------------------------------------! step 08

	summesion   = 0.0d0
	do i = 1,n
	summesion = summesion + t(i)*s(i)
	end do
	temp = summesion
	
	summesion   = 0.0d0
	do i = 1,n
	summesion = summesion + t(i)*t(i)
	end do
	
	omega = temp/summesion
	
!-----------------------------------------------------! step 09
	
	x = x + alpha*p + omega*s
	
!-----------------------------------------------------! step 11

	r = s - omega*t
	
!-----------------------------------------------------! step 10
	
	call norm
	
	it = it + 1
	
	end do !main loop
!*****************************************************!
	
	write(*,*) 'the solution is',x
	write(*,*) 'iteration required = ', it
	
end program BiCGSTAB


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

!-----------------------------------subroutine init---------------------------------------!

subroutine init
use var
implicit none
	integer :: i,j
	double precision :: summesion
	
	x = 0.0d0

!---------------------------------------------! step 01
	do i = 1,n
	summesion = 0.0d0
	do j = 1,n
	summesion = summesion + A(i,j)*x(j)
	end do
	r(i)  = b(i) - summesion
	end do
!---------------------------------------------! step 02	
	rs    = r
!---------------------------------------------! step 03
	rho   = 1.0d0
	alpha = 1.0d0
	omega = 1.0d0
!---------------------------------------------! step 04
	v = 0.0d0
	p = 0.0d0
!---------------------------------------------!

	call norm
	
end subroutine init

subroutine norm
use var
implicit none

	integer :: i
	double precision :: summesion
	!-------------------------------------------	
	
	summesion = 0.0d0
	
	do i = 1,n
	summesion = summesion + r(i)*r(i)
	end do
	
	norm_r = dsqrt(summesion)
	!-------------------------------------------
	summesion = 0.0d0
	
	do i = 1,n
	summesion = summesion + b(i)*b(i)
	end do
	
	norm_b = dsqrt(summesion)
end subroutine norm


!---------------------------------subroutine ab_def---------------------------------------!

subroutine ab_def
use var
implicit none

!---------------------first row

	A(1,1) = 1.0d0
	A(1,2) = 1.0d0
	A(1,3) = 1.0d0
	A(1,4) = 1.0d0

!---------------------end first row

	A(2,1) = 2.0d0
	A(2,2) = 3.0d0
	A(2,3) = 1.0d0
	A(2,4) = 1.0d0

	A(3,1) = 3.0d0
	A(3,2) = 4.0d0
	A(3,3) = 1.0d0
	A(3,4) = 1.0d0

	A(4,1) = 3.0d0
	A(4,2) = 4.0d0
	A(4,3) = 1.0d0
	A(4,4) = 2.0d0
	
	b(1) = 23.0d0
	b(2) = 13.0d0
	b(3) = 18.0d0
	b(4) = 23.0d0

	
end subroutine ab_def
My wiki