# CG, BICGSTAB(2) : problem with matrix operation and boundary conditions

 User Name Remember Me Password
 Register Blogs Members List Search Today's Posts Mark Forums Read

 February 12, 2010, 05:34 CG, BICGSTAB(2) : problem with matrix operation and boundary conditions #1 New Member     Join Date: May 2009 Posts: 9 Rep Power: 9 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``` In fact, for periodicity on x1, I copy point pl(Nx(1)-1) in ghost pl(0), and the point pl(2) in ghost pl(Nx(1)+1). 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, 19:47 #2 New Member     Join Date: May 2009 Posts: 9 Rep Power: 9 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``` GC : 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``` Type Grid : 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``` The CG do not converge, however the result is really near the solution : 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, 04:37 #3 New Member     Join Date: May 2009 Posts: 9 Rep Power: 9 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``` If used, you need to allocate x, and set it to 0.0 before the first time you use the routine, because the CG use the last solution calculated at the previous time as initial guess. This program is not very efficient, but it is simple and a may help to start with CG. With my best regards

 Thread Tools Display Modes Linear Mode

 Posting Rules You may not post new threads You may not post replies You may not post attachments You may not edit your posts BB code is On Smilies are On [IMG] code is On HTML code is OffTrackbacks are On Pingbacks are On Refbacks are On Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post Andy Chen FLUENT 4 July 10, 2009 01:20 pvc OpenFOAM Bugs 0 June 11, 2007 08:56

All times are GMT -4. The time now is 22:05.

 Contact Us - CFD Online - Top