[petsc-users] Using PETSC for a CFD solver
Mohammad Bahaa
m.bahaa.eldein at gmail.com
Tue Mar 18 07:44:05 CDT 2014
The "PETSC_COMM_SELF" worked like a charm for me, thanks alot
On Mon, Mar 10, 2014 at 4:50 PM, Matthew Knepley <knepley at gmail.com> wrote:
> 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
>
--
Mohamamd Bahaa ElDin
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20140318/15b1a028/attachment-0001.html>
More information about the petsc-users
mailing list