CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Main CFD Forum

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

Register Blogs Members List Search Today's Posts Mark Forums Read

Reply
 
LinkBack Thread Tools Display Modes
Old   February 12, 2010, 05:34
Default CG, BICGSTAB(2) : problem with matrix operation and boundary conditions
  #1
New Member
 
moomba's Avatar
 
Join Date: May 2009
Posts: 9
Rep Power: 8
moomba is on a distinguished road
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
moomba is offline   Reply With Quote

Old   February 12, 2010, 19:47
Default
  #2
New Member
 
moomba's Avatar
 
Join Date: May 2009
Posts: 9
Rep Power: 8
moomba is on a distinguished road
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
moomba is offline   Reply With Quote

Old   February 17, 2010, 04:37
Default
  #3
New Member
 
moomba's Avatar
 
Join Date: May 2009
Posts: 9
Rep Power: 8
moomba is on a distinguished road
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
moomba is offline   Reply With Quote

Reply

Thread Tools
Display Modes

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 Off
Trackbacks are On
Pingbacks are On
Refbacks are On


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


All times are GMT -4. The time now is 16:43.