 |
 |
 |
|
|
|
[Sponsors] | ||||
CG, BICGSTAB(2) : problem with matrix operation and boundary conditions |
 |
|
|
LinkBack | Thread Tools | Search this Thread | Display Modes |
February 12, 2010, 04:34
|
CG, BICGSTAB(2) : problem with matrix operation and boundary conditions
|
#1 |
|
New Member
Join Date: May 2009
Posts: 9
Rep Power: 16  |
Hi
I am trying to develop a Parallel Conjugate Gradient (CG) and a Parallel BICGSTAB(2) to solve Dynamic Pressure in a Low Mach DNS solver. I have been using the Van der Vorst paper.  The system is working well with perdiodic boundary conditions, but fails to work with Neumann (dP/dxn = 0, converge but bad solution) and Dirichlet (P=constant, do not converge). Because the problem is exactly the same for both system, I think it comes from the matrix multiplication. Could someone explain me where I made a mistake ? Here is my fortran 90 code for the BICGSTAB(2). The subroutine "Matrix_A_Multiply" is made here for periodicity on x1 and x2. The code is of course not optimised yet. Code:
module BICGSTAB_2
use typegrid
contains
subroutine BICGSTAB_2_CORE(Grid,b,x)
implicit none
type (grid_type) :: Grid
real(8), dimension(1:Grid%Nx(1),1:Grid%Nx(2),1:Grid%Nx(3)), intent(in) :: b
real(8), dimension(1:Grid%Nx(1),1:Grid%Nx(2),1:Grid%Nx(3)), intent(out) :: x
real(8), dimension(1:Grid%Nx(1),1:Grid%Nx(2),1:Grid%Nx(3)) :: xx0
real(8), dimension(1:Grid%Nx(1),1:Grid%Nx(2),1:Grid%Nx(3)) :: u,v,r,s,t,w
real(8) :: omega_1,omega_2,rho_1,rho_0,mu,nu,beta,gam,tau,alpha
integer :: i
xx0 = 0.0d0
CALL Matrix_A_Multiply(Grid,xx0,r) ! (x_in,x_out) --> x_out = A * x_in
r(:,:,:) = b(:,:,:) - r(:,:,:)
rho_0 = 1.0d0
u(:,:,:) = 0.0d0
alpha = 0.0d0
omega_2 = 1.0d0
do i=0,1000,2
rho_0 = - omega_2 * rho_0
! Even BiCG step :
rho_1 = Scalar_Product(Grid,r,r)
beta = alpha * rho_1 / rho_0
rho_0 = rho_1
u(:,:,:) = r(:,:,:) - beta * u(:,:,:)
CALL Matrix_A_Multiply(Grid,u,v)
gam = Scalar_Product(Grid,v,r) ! return v(1)*r(1) + v(2)*r(2) + v(3)*r(3) + ... + v(n)*r(n)
alpha = rho_0 / gam
r(:,:,:) = r(:,:,:) - alpha * v(:,:,:)
CALL Matrix_A_Multiply(Grid,r,s)
x(:,:,:) = x(:,:,:) + alpha * u(:,:,:)
! Odd BiCG step :
rho_1 = Scalar_Product(Grid,r,s)
beta = alpha * rho_1 / rho_0
rho_0 = rho_1
v(:,:,:) = s(:,:,:) - beta * v(:,:,:)
CALL Matrix_A_Multiply(Grid,v,w)
gam = Scalar_Product(Grid,w,r)
alpha = rho_0 / gam
u(:,:,:) = r(:,:,:) - beta * u(:,:,:)
r(:,:,:) = r(:,:,:) - alpha * v(:,:,:)
s(:,:,:) = s(:,:,:) - alpha * w(:,:,:)
CALL Matrix_A_Multiply(Grid,s,t)
! GCR(2)-part :
omega_1 = Scalar_Product(Grid,r,s)
mu = Scalar_Product(Grid,s,s)
nu = Scalar_Product(Grid,s,t)
tau = Scalar_Product(Grid,t,t)
omega_2 = Scalar_Product(Grid,r,t)
tau = tau - ( nu * nu / mu )
omega_2 = ( omega_2 - nu * omega_1 / mu ) / tau
! OPTI ?
omega_1 = ( omega_1 - nu * omega_2 ) / mu
x(:,:,:) = x(:,:,:) + omega_1 * r(:,:,:) + omega_2 * s(:,:,:) + alpha * u(:,:,:)
r(:,:,:) = r(:,:,:) - omega_1 * s(:,:,:) - omega_2 * t(:,:,:)
! test of convergence
if( max( ABS(maxval(r)) , ABS(minval(r)) ) < 1e-10 ) then
print *,"Converged in",i,"iterations"
return
end if
u(:,:,:) = u(:,:,:) - omega_1 * v(:,:,:) - omega_2 * w(:,:,:)
end do
print *,"ERROR, BICGSTAB(2) not converged"
end subroutine BICGSTAB_2_CORE
subroutine Matrix_A_Multiply(Grid,p,q)
implicit none
type (grid_type) :: Grid
real(8), dimension(1:Grid%Nx(1),1:Grid%Nx(2),1:Grid%Nx(3)), intent(in) :: p
real(8), dimension(1:Grid%Nx(1),1:Grid%Nx(2),1:Grid%Nx(3)), intent(out) :: q
real(8), dimension(0:Grid%Nx(1)+1,0:Grid%Nx(2)+1) :: pl
real(8), dimension(1:3) :: SDx
integer :: i,j,n
pl=77.0
pl(1:Grid%Nx(1),1:Grid%Nx(2)) = p(1:Grid%Nx(1),1:Grid%Nx(2),1)
! Opti to be made
pl( 0,1:Grid%Nx(2)) = pl(Grid%Nx(1)-1,1:Grid%Nx(2))
pl(Grid%Nx(1)+1,1:Grid%Nx(2)) = pl( 2,1:Grid%Nx(2))
pl(1:Grid%Nx(1), 0) = pl(1:Grid%Nx(1),Grid%Nx(2)-1)
pl(1:Grid%Nx(1),Grid%Nx(2)+1) = pl(1:Grid%Nx(1), 2)
do n=1,2
SDx(n)=(-1.0/(Grid%h(n)*Grid%h(n)))
end do
do j=1,Grid%Nx(2)
do i=1,Grid%Nx(1)
q(i,j,1)=SDx(1)*(-pl(i-1,j)-pl(i+1,j)+2.0*pl(i,j)) + &
& SDx(2)*(-pl(i,j-1)-pl(i,j+1)+2.0*pl(i,j))
end do
end do
return
end subroutine Matrix_A_Multiply
function Scalar_Product(Grid,x1,x2)
implicit none
type (grid_type) :: Grid
real(8), dimension(1:Grid%Nx(1),1:Grid%Nx(2),1:Grid%Nx(3)), intent(in) :: x1,x2
real(8) :: Scalar_Product
integer :: i,j,k
Scalar_Product=0.0d0
do k=1,Grid%Nx(3)
do j=1,Grid%Nx(2)
do i=1,Grid%Nx(1)
Scalar_Product = Scalar_Product + x1(i,j,k) * x2(i,j,k)
end do
end do
end do
return
end function Scalar_Product
end module BICGSTAB_2
For Neumann dP/dxn=0, I copy point pl(1) in ghost pl(0), and pl(Nx(1)) in ghost pl(Nx(1)+1). And for dirichlet P=0 on Nx(1) for example, I put pl(Nx(1))=0 and the output q(Nx(1))=0. I of course made a mistake, but where 
|
|
|
|
|
|
February 12, 2010, 18:47
|
|
#2 |
|
New Member
Join Date: May 2009
Posts: 9
Rep Power: 16 |
I made the problem far more simple : a GC to solve f=x² on [0,1] with Dirichlet on x (0 at x=0, and 1 at x=1
 ), and with dp/dy=0 on y border :Main : Code:
program Main
use typegrid
use CG
implicit none
type (grid_type) :: Grid
real(8), allocatable :: F(:,:),P(:,:)
real(8) :: PI,xx,yy
integer :: i,j
PI=3.141592653589793
Grid%Ndim=2
Grid%Nx(1)=100
Grid%Nx(2)=100
Grid%Nx(3)=1
Grid%h(1)=1.0/(Grid%Nx(1)-1)
Grid%h(2)=1.0/(Grid%Nx(2)-1)
Grid%h(3)=0.0
ALLOCATE(F(1:Grid%Nx(1),1:Grid%Nx(2)),P(1:Grid%Nx(1),1:Grid%Nx(2)))
F(:,:)=2.0
open(10,file="GC0.dat")
do i=1,Grid%Nx(1)
do j=1,Grid%Nx(2)
Write(10,*)(i-1)*Grid%h(1),(j-1)*Grid%h(2),F(i,j)
end do
end do
close(10)
CALL CG_CORE(Grid,F,P)
open(10,file="GC1.dat")
do i=1,Grid%Nx(1)
do j=1,Grid%Nx(2)
xx=(i-1)*Grid%h(1)
yy=(j-1)*Grid%h(2)
Write(10,*)(i-1)*Grid%h(1),(j-1)*Grid%h(2),P(i,j)
end do
end do
close(10)
DEALLOCATE(F,P)
end program Main
Code:
module CG
use typegrid
integer :: n01,n11,n02,n12
contains
subroutine CG_CORE(Grid,b,x)
implicit none
type (grid_type) :: Grid
real(8), dimension(1:Grid%Nx(1),1:Grid%Nx(2)), intent(in) :: b
real(8), dimension(1:Grid%Nx(1),1:Grid%Nx(2)), intent(out) :: x
real(8), dimension(1:Grid%Nx(1),1:Grid%Nx(2)) :: r,p,q,xx0
real(8) :: rho,rhop1,beta,alpha,res
integer :: i
n01 = 2
n11 = Grid%Nx(1)-1
n02 = 1
n12 = Grid%Nx(2)
xx0(n01:n11,n02:n12) = 0.0d0
CALL Matrix_A_Multiply(Grid,xx0,r) ! (x_in,x_out) --> x_out = A * x_in
r(n01:n11,n02:n12) = b(n01:n11,n02:n12) - r(n01:n11,n02:n12)
p(n01:n11,n02:n12) = 0.0d0
beta = 0.0d0
rho = Scalar_Product(Grid,r,r) ! return v(1)*r(1) + v(2)*r(2) + v(3)*r(3) + ... + v(n)*r(n)
do i=1,1000
p(n01:n11,n02:n12) = r(n01:n11,n02:n12) + beta * p(n01:n11,n02:n12)
CALL Matrix_A_Multiply(Grid,p,q)
alpha = rho / Scalar_Product(Grid,p,q)
x(n01:n11,n02:n12) = x(n01:n11,n02:n12) + alpha * p(n01:n11,n02:n12)
r(n01:n11,n02:n12) = r(n01:n11,n02:n12) - alpha * q(n01:n11,n02:n12)
res = Residut(Grid,r)
if( res < 1e-6 ) then
print *, "GC converged in",i,"IT with a residut of",res
return
end if
! print *,res
rhop1 = Scalar_Product(Grid,r,r)
beta = rhop1 / rho
rho = rhop1
end do
print *,"ERROR, CG not converged"
end subroutine CG_CORE
subroutine Matrix_A_Multiply(Grid,p,q)
implicit none
type (grid_type) :: Grid
real(8), dimension(1:Grid%Nx(1),1:Grid%Nx(2)), intent(in) :: p
real(8), dimension(1:Grid%Nx(1),1:Grid%Nx(2)), intent(out) :: q
real(8), dimension(0:Grid%Nx(1)+1,0:Grid%Nx(2)+1) :: pl ! copy with ghost, not optimized
real(8), dimension(1:3) :: SDx
integer :: i,j,n
pl=77.0 ! dummy value
pl(1:Grid%Nx(1),1:Grid%Nx(2)) = p(1:Grid%Nx(1),1:Grid%Nx(2))
! Dirichlet : 0 en x=0, et 1 en x=1
pl( 1,1:Grid%Nx(2)) = 0.0!(( 1-1)*Grid%h(1))**2
pl(Grid%Nx(1),1:Grid%Nx(2)) = 1.0!((Grid%Nx(1)-1)*Grid%h(1))**2
! Neumann dp/dn = 0
pl(1:Grid%Nx(1), 0) = pl(1:Grid%Nx(1), 1)
pl(1:Grid%Nx(1),Grid%Nx(2)+1) = pl(1:Grid%Nx(1),Grid%Nx(2))
do n=1,2
SDx(n)=(1.0d0/(Grid%h(n)*Grid%h(n)))
end do
do j=n02,n12
do i=n01,n11
q(i,j) = SDx(1)*(pl(i-1,j)+pl(i+1,j)-2.0d0*pl(i,j)) + &
& SDx(2)*(pl(i,j-1)+pl(i,j+1)-2.0d0*pl(i,j))
end do
end do
return
end subroutine Matrix_A_Multiply
function Residut(Grid,r)
implicit none
type (grid_type) :: Grid
real(8), dimension(1:Grid%Nx(1),1:Grid%Nx(2)), intent(in) :: r
real(8) :: Residut
integer :: i,j
Residut=0.0d0
do j=n02,n12
do i=n01,n11
Residut = max ( Residut, ABS(r(i,j)) )
end do
end do
return
end function Residut
function Scalar_Product(Grid,x1,x2)
implicit none
type (grid_type) :: Grid
real(8), dimension(1:Grid%Nx(1),1:Grid%Nx(2)), intent(in) :: x1,x2
real(8) :: Scalar_Product
integer :: i,j
Scalar_Product=0.0d0
do j=n02,n12
do i=n01,n11
Scalar_Product = Scalar_Product + x1(i,j) * x2(i,j)
end do
end do
return
end function Scalar_Product
end module CG
Code:
module typegrid
implicit none
type grid_type
integer :: Ndim
integer :: NPL
integer, dimension(1:3) :: Nx
real, dimension(1:3) :: X0
real, dimension(1:3) :: Lx
real, dimension(1:3) :: h
real, dimension(1:100,1:3) :: XP
end type grid_type
end module typegrid
 As you can see, the value is a little different near x=1... Does anyone have an idée ? It would really help me
|
|
|
|
|
|
|
February 17, 2010, 03:37
|
|
#3 |
|
New Member
Join Date: May 2009
Posts: 9
Rep Power: 16 |
Problem solve, here are the sources if it can help someone.
It is a simple conjugate gradient with neumann, periodic and dirichlet boundary conditions. Code:
! Neumann at face 1.1 BoundCG(1,1) = 1 ! Dirichlet at face 1.2 BoundCG(1,2) = -1 PRES(Grid%Nx(1),:,:) = 0.0d0 ! periodic on axe y BoundCG(2,1) = 12 BoundCG(2,2) = 12 CALL CG_CORE(Grid,PRES,BoundCG) Code:
module CG
use typegrid
contains
subroutine CG_CORE(Grid,b,BoundCG)
use commonGC
implicit none
type (grid_type) :: Grid
real(8), dimension(1:Grid%Nx(1),1:Grid%Nx(2),1:Grid%Nx(3)), intent(inout) :: b
integer, dimension(1:3,1:2), intent(in) :: BoundCG
real(8), dimension(1:Grid%Nx(1),1:Grid%Nx(2),1:Grid%Nx(3)) :: r,p,q,w
real(8), dimension(1:3) :: SDx
real(8) :: rho,rhop1,beta,alpha,res,a
integer :: i,n
! BoundCG
! 10 = simple com
! 12 = periodic com
! 1 = inflow (dP/dn = 0)
! -1 = outflow (P = fixed)
! Dirichlet Boundary Conditions, set fixed value in initial guess, this values will not be overwritten during the CG process
if(BoundCG(1,1) == -1) x(1,:,:) = b(1,:,:)
if(BoundCG(1,2) == -1) x(Grid%Nx(1),:,:) = b(Grid%Nx(1),:,:)
if(Grid%ndim > 1) then
if(BoundCG(2,1) == -1) x(:,1,:) = b(:,1,:)
if(BoundCG(2,2) == -1) x(:,Grid%Nx(2),:) = b(:,Grid%Nx(2),:)
end if
if(Grid%ndim > 2) then
if(BoundCG(3,1) == -1) x(:,:,1) = b(:,:,1)
if(BoundCG(3,2) == -1) x(:,:,Grid%Nx(3)) = b(:,:,Grid%Nx(3))
end if
CALL Matrix_A_Multiply(Grid,x,r,BoundCG) ! (x_in,x_out) --> x_out = A * x_in
r(:,:,:) = b(:,:,:) - r(:,:,:)
! Dirichlet Boundary Conditions, set residut at 0.0 (value is already correct)
if(BoundCG(1,1) == -1) r(1,:,:) = 0.0d0
if(BoundCG(1,2) == -1) r(Grid%Nx(1),:,:) = 0.0d0
if(Grid%ndim > 1) then
if(BoundCG(2,1) == -1) r(:,1,:) = 0.0d0
if(BoundCG(2,2) == -1) r(:,Grid%Nx(2),:) = 0.0d0
end if
if(Grid%ndim > 2) then
if(BoundCG(3,1) == -1) r(:,:,1) = 0.0d0
if(BoundCG(3,2) == -1) r(:,:,Grid%Nx(3)) = 0.0d0
end if
p(:,:,:) = 0.0d0
beta = 0.0d0
! PRECONDITIONNEMENT
do n=1,Grid%ndim
SDx(n)=(1.0d0/(Grid%h(n)*Grid%h(n)))
end do
a=1.0d0/(-2.0d0*SDx(1)-2.0d0*SDx(2))
w(:,:,:) = a * r(:,:,:)
rho = Scalar_Product(Grid,r,w) ! return v(1)*r(1) + v(2)*r(2) + v(3)*r(3) + ... + v(n)*r(n)
do i=1,2000
p(:,:,:) = w(:,:,:) + beta * p(:,:,:)
CALL Matrix_A_Multiply(Grid,p,q,BoundCG)
alpha = rho / Scalar_Product(Grid,p,q)
x(:,:,:) = x(:,:,:) + alpha * p(:,:,:)
r(:,:,:) = r(:,:,:) - alpha * q(:,:,:)
res = Residut(Grid,r)
if( res < 1e-6 ) then
print *, "GC converged in",i,"IT with a residut of",res
b(:,:,:) = x(:,:,:)
return
end if
! print *,res
w(:,:,:) = a * r(:,:,:)
rhop1 = Scalar_Product(Grid,r,w)
beta = rhop1 / rho
rho = rhop1
end do
print *,"ERROR, CG not converged"
end subroutine CG_CORE
subroutine Matrix_A_Multiply(Grid,p,q,BoundCG)
implicit none
type (grid_type) :: Grid
real(8), dimension(1:Grid%Nx(1),1:Grid%Nx(2),1:Grid%Nx(3)), intent(in) :: p
real(8), dimension(1:Grid%Nx(1),1:Grid%Nx(2),1:Grid%Nx(3)), intent(out) :: q
integer, dimension(1:3,1:2), intent(in) :: BoundCG
real(8), dimension(0:Grid%Nx(1)+1,0:Grid%Nx(2)+1,0:Grid%Nx(3)+1) :: pl ! copy with ghost, not optimized
real(8), dimension(1:3) :: SDx
integer :: i,j,n
pl=77.0
pl(1:Grid%Nx(1),1:Grid%Nx(2),1) = p(1:Grid%Nx(1),1:Grid%Nx(2),1)
! Neumann Boundaries, dP/dn = 0 , not optimized (should not copy, but simply not compute)
if(BoundCG(1,1) == 1) pl(0,:,:) = pl(2,:,:)
if(BoundCG(1,2) == 1) pl(Grid%Nx(1)+1,:,:) = pl(Grid%Nx(1)-1,:,:)
if(Grid%ndim > 1) then
if(BoundCG(2,1) == 1) pl(:,0,:) = pl(:,2,:)
if(BoundCG(2,2) == 1) pl(:,Grid%Nx(2)+1,:) = pl(:,Grid%Nx(2)-1,:)
end if
if(Grid%ndim > 2) then
if(BoundCG(3,1) == 1) pl(:,:,0) = pl(:,:,2)
if(BoundCG(3,2) == 1) pl(:,:,Grid%Nx(3)+1) = pl(:,:,Grid%Nx(3)-1)
end if
! Periodic Boundaries
if(BoundCG(1,1) == 12) pl(0,:,:) = pl(Grid%Nx(1)-1,:,:)
if(BoundCG(1,2) == 12) pl(Grid%Nx(1)+1,:,:) = pl(2,:,:)
if(Grid%ndim > 1) then
if(BoundCG(2,1) == 12) pl(:,0,:) = pl(:,Grid%Nx(2)-1,:)
if(BoundCG(2,2) == 12) pl(:,Grid%Nx(2)+1,:) = pl(:,2,:)
end if
if(Grid%ndim > 2) then
if(BoundCG(3,1) == 12) pl(:,:,0) = pl(:,:,Grid%Nx(3)-1)
if(BoundCG(3,2) == 12) pl(:,:,Grid%Nx(3)+1) = pl(:,:,2)
end if
do n=1,2
SDx(n)=(1.0/(Grid%h(n)*Grid%h(n)))
end do
do j=1,Grid%Nx(2)
do i=1,Grid%Nx(1)
q(i,j,1)=SDx(1)*(pl(i-1,j,1)+pl(i+1,j,1)-2.0*pl(i,j,1)) + &
& SDx(2)*(pl(i,j-1,1)+pl(i,j+1,1)-2.0*pl(i,j,1))
end do
end do
! Dirichelt boundaries, set q at 0 (to not alter the x fixed value)
if(BoundCG(1,1) == -1) q(1,:,:) = 0.0d0
if(BoundCG(1,2) == -1) q(Grid%Nx(1),:,:) = 0.0d0
if(Grid%ndim > 1) then
if(BoundCG(2,1) == -1) q(:,1,:) = 0.0d0
if(BoundCG(2,2) == -1) q(:,Grid%Nx(2),:) = 0.0d0
end if
if(Grid%ndim > 2) then
if(BoundCG(3,1) == -1) q(:,:,1) = 0.0d0
if(BoundCG(3,2) == -1) q(:,:,Grid%Nx(3)) = 0.0d0
end if
return
end subroutine Matrix_A_Multiply
function Residut(Grid,r)
implicit none
type (grid_type) :: Grid
real(8), dimension(1:Grid%Nx(1),1:Grid%Nx(2),1:Grid%Nx(3)), intent(in) :: r
real(8) :: Residut
integer :: i,j,k
Residut=0.0d0
do k=1,Grid%Nx(3)
do j=1,Grid%Nx(2)
do i=1,Grid%Nx(1)
Residut = max ( Residut, ABS(r(i,j,k)) )
end do
end do
end do
return
end function Residut
function Scalar_Product(Grid,x1,x2)
implicit none
type (grid_type) :: Grid
real(8), dimension(1:Grid%Nx(1),1:Grid%Nx(2),1:Grid%Nx(3)), intent(in) :: x1,x2
real(8) :: Scalar_Product
integer :: i,j,k
Scalar_Product=0.0d0
do k=1,Grid%Nx(3)
do j=1,Grid%Nx(2)
do i=1,Grid%Nx(1)
Scalar_Product = Scalar_Product + x1(i,j,k) * x2(i,j,k)
end do
end do
end do
return
end function Scalar_Product
end module CG
Code:
module commonGC implicit none real(8), allocatable, dimension(:,:,:), save :: x end module commonGC This program is not very efficient, but it is simple and a may help to start with CG. With my best regards |
|
|
|
|
|
|
|
|
Similar Threads
|
||||
| Thread | Thread Starter | Forum | Replies | Last Post |
| open channel problem (free surface) | Andy Chen | FLUENT | 4 | July 10, 2009 01:20 |
| Problem with boundary file | pvc | OpenFOAM Bugs | 0 | June 11, 2007 08:56 |
Linear Mode
