[petsc-users] Using PETSC for a CFD solver
Mohammad Bahaa
m.bahaa.eldein at gmail.com
Mon Mar 10 09:36:47 CDT 2014
Hi,
I'm pretty new to PETSC, so pardon me if the question is primitive somehow,
I used *METIS *to partition my grid (represented by a system of linear
equations Ax=b) to a number of sub-grids, say 4 sub-grids, with 4 different
systems of linear Equations (A1x1=b1, A2x2=b2, ...), can anyone post an
example showing how to solve these "n" sub-systems simultaneously, I've
tried the following program, but it's not working correctly, as when I
use *MatGetOwnershipRange
*in each process I find that A1 ownership range is 1/4 the matrix size for
the first process, while it should be all of it.
subroutine test_drive_2
implicit none
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! Include files
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
! This program uses CPP for preprocessing, as indicated by the use of
! PETSc include files in the directory petsc/include/finclude. This
! convention enables use of the CPP preprocessor, which allows the use
! of the #include statements that define PETSc objects and variables.
!
! Use of the conventional Fortran include statements is also supported
! In this case, the PETsc include files are located in the directory
! petsc/include/foldinclude.
!
! Since one must be very careful to include each file no more than once
! in a Fortran routine, application programmers must exlicitly list
! each file needed for the various PETSc components within their
! program (unlike the C/C++ interface).
!
! See the Fortran section of the PETSc users manual for details.
!
! The following include statements are required for KSP Fortran programs:
! petscsys.h - base PETSc routines
! petscvec.h - vectors
! petscmat.h - matrices
! petscksp.h - Krylov subspace methods
! petscpc.h - preconditioners
! Other include statements may be needed if using additional PETSc
! routines in a Fortran program, e.g.,
! petscviewer.h - viewers
! petscis.h - index sets
!
! main includes
#include <finclude/petscsys.h>
#include <finclude/petscvec.h>
#include <finclude/petscmat.h>
#include <finclude/petscksp.h>
#include <finclude/petscpc.h>
! includes for F90 specific functions
#include <finclude/petscvec.h90>
#include <finclude/petscmat.h90>
#include <finclude/petscksp.h90>
#include <finclude/petscpc.h90>
!
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! Variable declarations
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
! Variables:
! ksp - linear solver context
! ksp - Krylov subspace method context
! pc - preconditioner context
! x, b, u - approx solution, right-hand-side, exact solution vectors
! A - matrix that defines linear system
! its - iterations for convergence
! norm - norm of error in solution
!
Vec x,b,u
Mat A
KSP ksp
PC pc
PetscReal norm,tol
PetscErrorCode ierr
PetscInt i,n,col(3),its,i1,i2,i3
PetscBool flg
PetscMPIInt size,rank
PetscScalar none,one,value(3), testa
PetscScalar, pointer :: xx_v(:)
PetscScalar, allocatable, dimension(:) :: myx
PetscOffset i_x
!real(4) :: myx(10), myu(10), myb(10)
!real(8), allocatable, dimension(:) :: myx
integer :: ic, nc, ncmax, nz, ncols, j
integer :: fileunit, ione
integer, allocatable, dimension(:,:) :: neighb, cols
integer, allocatable, dimension(:) :: nnz, vcols
real(8), allocatable, dimension(:,:) :: acoef, vals
real(8), allocatable, dimension(:) :: ap, su, vvals
character (len=100) :: rankstring, filename, folder
real(8) :: atol, rtol, dtol
integer :: mxit, istart, iend
real(8) :: rvar, minvalx
call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)
call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! Load the linear system
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!mbs read data file for experimentation
!if(rank.EQ.0) read(*,'(A)'), folder
folder = '60'
write(rankstring,'(I)'), rank
filename = trim(adjustl(folder)) // '/linearsys_' //
trim(adjustl(rankstring)) // '.txt'
fileunit = 9000 + rank
open(unit=fileunit, file=trim(filename))
read(fileunit,*), ncmax
read(fileunit,*), nc
read(fileunit,*),
allocate( neighb(6,ncmax), acoef(6,ncmax), ap(ncmax), su(ncmax) )
allocate( nnz(0:ncmax-1), cols(6,0:ncmax-1), vals(6,0:ncmax-1) )
allocate( vcols(0:ncmax-1), vvals(0:ncmax-1) )
!allocate( xx_v(ncmax), myx(ncmax) )
!allocate( xx_v(0:ncmax), myx(0:ncmax) )
allocate( myx(ncmax) )
do i=1,nc
read(fileunit,'(I)'), ic
read(fileunit,'(6I)'), ( neighb(j,ic), j=1,6 )
read(fileunit,'(6F)'), ( acoef(j,ic), j=1,6 )
read(fileunit,'(2F)'), ap(ic), su(ic)
read(fileunit,*),
enddo
close(fileunit)
nz = 7
nnz = 0
do ic=0,nc-1
! values for coefficient matrix (diagonal)
nnz(ic) = nnz(ic) + 1
cols(nnz(ic),ic) = ic
vals(nnz(ic),ic) = ap(ic+1)
! values for coefficient matrix (off diagonal)
do j=1,6
if(neighb(j,ic+1).GT.0)then
nnz(ic) = nnz(ic) + 1
cols(nnz(ic),ic) = neighb(j,ic+1) - 1
vals(nnz(ic),ic) = acoef(j,ic+1)
endif
enddo
! values for RHS
vcols(ic) = ic
vvals(ic) = su(ic+1)
enddo
! add dummy values for remaining rows (if any)
if(ncmax.GT.nc)then
do ic=nc,ncmax-1
! coeff matrix
nnz(ic) = 1
cols(nnz(ic),ic) = ic
vals(nnz(ic),ic) = 1.0
! RHS
vcols(ic) = ic
vvals(ic) = 0.0
enddo
endif
print*, 'rank', rank, 'says nc is', nc
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! Beginning of program
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!if (size .ne. 1) then
! call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
! if (rank .eq. 0) then
! write(6,*) 'This is a uniprocessor example only!'
! endif
! SETERRQ(PETSC_COMM_WORLD,1,' ',ierr)
!endif
ione = 1
none = -1.0
one = 1.0
n = ncmax
i1 = 1
i2 = 2
i3 = 3
call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-n',n,flg,ierr)
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! Compute the matrix and right-hand-side vector that define
! the linear system, Ax = b.
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! Create matrix. When using MatCreate(), the matrix format can
! be specified at runtime.
call MatCreate(PETSC_COMM_WORLD,A,ierr)
!call MatCreateSeqAij(PETSC_COMM_SELF,A,ierr)
call MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n,ierr)
call MatSetFromOptions(A,ierr)
call MatSetUp(A,ierr)
call MatGetOwnershipRange(A,istart,iend,ierr)
print*, rank, istart, iend
! Assemble matrix.
! - Note that MatSetValues() uses 0-based row and column numbers
! in Fortran as well as in C (as set here in the array "col").
! value(1) = -1.0
! value(2) = 2.0
! value(3) = -1.0
! do 50 i=1,n-2
! col(1) = i-1
! col(2) = i
! col(3) = i+1
! call MatSetValues(A,i1,i,i3,col,value,INSERT_VALUES,ierr)
!50 continue
! i = n - 1
! col(1) = n - 2
! col(2) = n - 1
! call MatSetValues(A,i1,i,i2,col,value,INSERT_VALUES,ierr)
! i = 0
! col(1) = 0
! col(2) = 1
! value(1) = 2.0
! value(2) = -1.0
! call MatSetValues(A,i1,i,i2,col,value,INSERT_VALUES,ierr)
! call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr)
! call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr)
do ic=0,ncmax-1
call
MatSetValues(A,ione,ic,nnz(ic),cols(1:nnz(ic),ic),vals(1:nnz(ic),ic),INSERT_VALUES,ierr)
enddo
call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr)
call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr)
! Create vectors. Note that we form 1 vector from scratch and
! then duplicate as needed.
call VecCreate(PETSC_COMM_WORLD,x,ierr)
!call VecCreateSeq(PETSC_COMM_SELF,x,ierr)
call VecSetSizes(x,PETSC_DECIDE,n,ierr)
call VecSetFromOptions(x,ierr)
call VecDuplicate(x,b,ierr)
call VecDuplicate(x,u,ierr)
! Set exact solution; then compute right-hand-side vector.
call VecSet(u,one,ierr)
call VecSet(x,one*10,ierr)
!call MatMult(A,u,b,ierr)
! set source terms vector
call VecSetValues(b,ncmax,vcols,vvals,INSERT_VALUES,ierr)
call VecAssemblyBegin(b,ierr)
call VecAssemblyEnd(b,ierr)
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! Create the linear solver and set various options
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! Create linear solver context
call KSPCreate(PETSC_COMM_WORLD,ksp,ierr)
! Set operators. Here the matrix that defines the linear system
! also serves as the preconditioning matrix.
call KSPSetOperators(ksp,A,A,DIFFERENT_NONZERO_PATTERN,ierr)
! Set linear solver defaults for this problem (optional).
! - By extracting the KSP and PC contexts from the KSP context,
! we can then directly directly call any KSP and PC routines
! to set various options.
! - The following four statements are optional; all of these
! parameters could alternatively be specified at runtime via
! KSPSetFromOptions();
call KSPGetPC(ksp,pc,ierr)
call PCSetType(pc,PCJACOBI,ierr)
atol = 1.d-12
rtol = 1.d-12
dtol = 1.d10
mxit = 100
! call KSPSetTolerances(ksp,tol,PETSC_DEFAULT_DOUBLE_PRECISION, &
! & PETSC_DEFAULT_DOUBLE_PRECISION,PETSC_DEFAULT_INTEGER,ierr)
call KSPSetTolerances(ksp,atol,rtol,dtol,mxit,ierr)
! Set runtime options, e.g.,
! -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
! These options will override those specified above as long as
! KSPSetFromOptions() is called _after_ any other customization
! routines.
call KSPSetFromOptions(ksp,ierr)
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! Solve the linear system
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
call KSPSetType(ksp,KSPBCGS,ierr)
call KSPSolve(ksp,b,x,ierr)
! View solver info; we could instead use the option -ksp_view
!call KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD,ierr)
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! Check solution and clean up
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
call VecGetArrayF90(x,xx_v,ierr)
!xx_v = 5.1d0
!do ic=1,ncmax
! myx(ic) = xx_v(ic)
!enddo
myx = xx_v
!call VecGetArray(x,myx,i_x,ierr)
!value = x_array(i_x + 1)
!call VecRestoreArray(x,myx,i_x,ierr)
!rvar = xx_v(3)
call VecRestoreArrayF90(x,xx_v,ierr)
!
!call VecGetArrayF90(b,xx_v,ierr)
!myb = xx_v
!call VecRestoreArrayF90(x,xx_v,ierr)
!
!call VecView(x,PETSC_VIEWER_STDOUT_SELF)
!call MatView(a,PETSC_VIEWER_STDOUT_SELF)
!print*, 'rank', rank, 'says max x is', maxval(myx)
! print*, xx_v
! Check the error
call MatMult(A,x,u,ierr)
call VecAXPY(u,none,b,ierr)
call VecNorm(u,NORM_2,norm,ierr)
call KSPGetIterationNumber(ksp,its,ierr)
if (norm .gt. 1.e-12) then
write(6,100) norm,its
else
write(6,200) its
endif
100 format('Norm of error = ',e11.4,', Iterations = ',i5)
200 format('Norm of error < 1.e-12,Iterations = ',i5)
!call KSPGetSolution(ksp,myx)
minvalx = 1.0e15
do ic=1,ncmax
if(myx(ic).LT.minvalx) minvalx = myx(ic)
enddo
!write(*,300), rank, maxval(myx(1:nc)), minvalx
if(rank.EQ.0) print*, myx
300 format('Rank ', I1, ' says max/min x are: ', F12.4, ' / ', F12.4)
! Free work space. All PETSc objects should be destroyed when they
! are no longer needed.
call VecDestroy(x,ierr)
call VecDestroy(u,ierr)
call VecDestroy(b,ierr)
call MatDestroy(A,ierr)
call KSPDestroy(ksp,ierr)
call PetscFinalize(ierr)
continue
end
--
Mohamamd Bahaa ElDin
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20140310/6f11b8d3/attachment-0001.html>
More information about the petsc-users
mailing list