[moomba]

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.

      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