[petsc-users] Using PETSC for a CFD solver
Matthew Knepley
knepley at gmail.com
Mon Mar 10 09:50:19 CDT 2014
On Mon, Mar 10, 2014 at 9:36 AM, Mohammad Bahaa <m.bahaa.eldein at gmail.com>wrote:
> 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.
>
I will answer this two ways. First, here is the "PETSc strategy" for doing
the same thing.
1) Write a code that assembles and solves the entire thing
2) Use -pc_type bjacobi -ksp_type preonly -sub_ksp_type <whatever you
want> -mat
3) You can use ParMetis inside PETSc with the MatPartitioning
This will solve the individual systems with no coupling (this is what it
sounds like you want above).
If you want to manage everything yourself, and you want to form individual
systems on every process,
just create the solvers using a smaller communicator. PETSC_COMM_SELF means
that every system
is serial. You can make smaller comms with MPI_Comm_split() if you want
smaller comms, but some
parallelism for each system.
Thanks,
Matt
> 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
>
--
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20140310/e4eaf9a9/attachment-0001.html>
More information about the petsc-users
mailing list