CFD Online Discussion Forums (http://www.cfd-online.com/Forums/)
-   Main CFD Forum (http://www.cfd-online.com/Forums/main/)
-   -   CG, BICGSTAB(2) : problem with matrix operation and boundary conditions (http://www.cfd-online.com/Forums/main/72622-cg-bicgstab-2-problem-matrix-operation-boundary-conditions.html)

 moomba February 12, 2010 05:34

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

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 :confused:

 moomba February 12, 2010 19:47

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 :D), 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 :

http://img534.imageshack.us/img534/4...twindowid0.png

As you can see, the value is a little different near x=1...

Does anyone have an idée ? It would really help me :)

 moomba February 17, 2010 04:37

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

 All times are GMT -4. The time now is 01:41.