TDMA.f90 - Solution of system of linear equations by Thomas method
From CFD-Wiki
(Difference between revisions)
(New page: <pre> !Sample program for solving Lid-driven cavity flow test using SIMPLE-algorithm ! solution of linear system of equations by Thomas algorithm modul !Copyright (C) 2010 Michail Kirič...) |
|||
Line 4: | Line 4: | ||
! solution of linear system of equations by Thomas algorithm modul | ! solution of linear system of equations by Thomas algorithm modul | ||
!Copyright (C) 2010 Michail Kiričkov | !Copyright (C) 2010 Michail Kiričkov | ||
+ | !Copyright (C) 2016 Michail Kiričkov, Kaunas University for Technology | ||
!This program is free software; you can redistribute it and/or | !This program is free software; you can redistribute it and/or | ||
Line 21: | Line 22: | ||
!********************************************************************** | !********************************************************************** | ||
Subroutine TDMA_1(NF) | Subroutine TDMA_1(NF) | ||
- | |||
include 'icomm_1.f90' | include 'icomm_1.f90' | ||
- | |||
DOUBLE PRECISION P(nx),Q(nx) | DOUBLE PRECISION P(nx),Q(nx) | ||
- | + | !----------------------------------------------------------- | |
- | ! | + | |
Do 101 J = 2, NYmax | Do 101 J = 2, NYmax | ||
- | |||
P(1) = 0. | P(1) = 0. | ||
- | + | Q(1) = F(1,j,nf) | |
- | + | P(NXmaxC) = 0. | |
- | + | Q(NXmaxC) = F(NXmaxC,j,nf) | |
- | + | ||
- | + | ||
! Forward Elimination | ! Forward Elimination | ||
- | + | Do 10 i = 2,NXmaxC-1 | |
- | + | ||
- | + | ||
temp = Ap(i,j,nf) - Aw(i,j) * P(i-1) | temp = Ap(i,j,nf) - Aw(i,j) * P(i-1) | ||
- | + | Spp= Sp(i,j,nf) + As(i,j) * F(i,j-1,nf) + & | |
- | + | An(i,j) * F(i,j+1,nf) | |
- | + | P(i) = Ae(i,j) / temp | |
- | + | Q(i) = (Spp + Aw(i,j)*Q(i-1)) / temp | |
- | + | 10 continue | |
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
! Back Substitution | ! Back Substitution | ||
- | |||
Do 20 i = NXmaxC-1,1,-1 | Do 20 i = NXmaxC-1,1,-1 | ||
- | |||
F(i,j,nf) = P(i)*F(i+1,j,nf) + Q(i) | F(i,j,nf) = P(i)*F(i+1,j,nf) + Q(i) | ||
- | |||
20 continue | 20 continue | ||
- | |||
101 continue | 101 continue | ||
- | ! | + | !----------------------------------------------------------- |
- | + | ||
Do 301 J = NYmax,2,-1 | Do 301 J = NYmax,2,-1 | ||
- | |||
P(1) = 0. | P(1) = 0. | ||
- | + | Q(1) = F(1,j,nf) | |
- | + | P(NXmaxC) = 0. | |
- | + | Q(NXmaxC) = F(NXmaxC,j,nf) | |
- | + | ||
- | + | ||
! Forward Elimination | ! Forward Elimination | ||
- | |||
Do 32 i = 2,NXmaxC-1 | Do 32 i = 2,NXmaxC-1 | ||
- | + | temp = Ap(i,j,nf) - Aw(i,j) * P(i-1) | |
- | + | Spp= Sp(i,j,nf) + As(i,j) * F(i,j-1,nf) + & | |
- | + | An(i,j) * F(i,j+1,nf) | |
- | + | P(i) = Ae(i,j) / temp | |
- | + | Q(i) = (Spp + Aw(i,j)*Q(i-1)) / temp | |
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
32 continue | 32 continue | ||
- | |||
! Back Substitution | ! Back Substitution | ||
- | |||
Do 30 i = NXmaxC-1,2,-1 | Do 30 i = NXmaxC-1,2,-1 | ||
- | |||
F(i,j,nf) = P(i)*F(i+1,j,nf) + Q(i) | F(i,j,nf) = P(i)*F(i+1,j,nf) + Q(i) | ||
- | |||
30 continue | 30 continue | ||
- | |||
301 continue | 301 continue | ||
- | !-------------------------------------------------------------- | + | !-------------------------------------------------------------- |
- | + | ||
- | + | ||
Return | Return | ||
End | End | ||
</pre> | </pre> |
Latest revision as of 15:00, 19 May 2016
!Sample program for solving Lid-driven cavity flow test using SIMPLE-algorithm ! solution of linear system of equations by Thomas algorithm modul !Copyright (C) 2010 Michail Kiričkov !Copyright (C) 2016 Michail Kiričkov, Kaunas University for Technology !This program is free software; you can redistribute it and/or !modify it under the terms of the GNU General Public License !as published by the Free Software Foundation; either version 2 !of the License, or (at your option) any later version. !This program is distributed in the hope that it will be useful, !but WITHOUT ANY WARRANTY; without even the implied warranty of !MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the !GNU General Public License for more details. !You should have received a copy of the GNU General Public License !along with this program; if not, write to the Free Software !Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. !********************************************************************** Subroutine TDMA_1(NF) include 'icomm_1.f90' DOUBLE PRECISION P(nx),Q(nx) !----------------------------------------------------------- Do 101 J = 2, NYmax P(1) = 0. Q(1) = F(1,j,nf) P(NXmaxC) = 0. Q(NXmaxC) = F(NXmaxC,j,nf) ! Forward Elimination Do 10 i = 2,NXmaxC-1 temp = Ap(i,j,nf) - Aw(i,j) * P(i-1) Spp= Sp(i,j,nf) + As(i,j) * F(i,j-1,nf) + & An(i,j) * F(i,j+1,nf) P(i) = Ae(i,j) / temp Q(i) = (Spp + Aw(i,j)*Q(i-1)) / temp 10 continue ! Back Substitution Do 20 i = NXmaxC-1,1,-1 F(i,j,nf) = P(i)*F(i+1,j,nf) + Q(i) 20 continue 101 continue !----------------------------------------------------------- Do 301 J = NYmax,2,-1 P(1) = 0. Q(1) = F(1,j,nf) P(NXmaxC) = 0. Q(NXmaxC) = F(NXmaxC,j,nf) ! Forward Elimination Do 32 i = 2,NXmaxC-1 temp = Ap(i,j,nf) - Aw(i,j) * P(i-1) Spp= Sp(i,j,nf) + As(i,j) * F(i,j-1,nf) + & An(i,j) * F(i,j+1,nf) P(i) = Ae(i,j) / temp Q(i) = (Spp + Aw(i,j)*Q(i-1)) / temp 32 continue ! Back Substitution Do 30 i = NXmaxC-1,2,-1 F(i,j,nf) = P(i)*F(i+1,j,nf) + Q(i) 30 continue 301 continue !-------------------------------------------------------------- Return End