Sample code for BiCGSTAB - Fortran 90
From CFD-Wiki
(Difference between revisions)
Line 3: | Line 3: | ||
! ! | ! ! | ||
! Copyright (c) 2015, Biswajit Ghosh ! | ! Copyright (c) 2015, Biswajit Ghosh ! | ||
- | ! mebiswajithit@gmail.com | + | ! mebiswajithit@gmail.com ! |
! All rights reserved. ! | ! All rights reserved. ! | ||
! ! | ! ! |
Revision as of 16:33, 13 September 2016
!-----------------------------------MIT LICENCE-----------------------------------------! ! ! ! Copyright (c) 2015, Biswajit Ghosh ! ! mebiswajithit@gmail.com ! ! 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