On Tue, Feb 9, 2016 at 9:10 AM, Hom Nath Gharti <hng.em...@gmail.com> wrote:
> Thank you so much Barry! > > For my small test case, with -pc_fieldsplit_block_size 4, the program > runs, although the answer was not correct. At least now I get > something to look upon. I am using PCFieldSplitSetIS to set the > fields. Do I still need to use -pc_fieldsplit_block_size? > No, this is only for co-located discretizations. > In my case each grid point may have different variable sets. For > example, the grid point in the solid region has displacement and > gravity potential: ux, uy, uz, and \phi. Whereas the grid point in the > fluid region has scalar potential, pressure and gravity potential: > \chi, p, and \phi. And the solid-fluid interface has all of them. Do > you think I can still use PCFIELDSPLIT in this situation? > We have examples, like SNES ex62, which do exactly this. Thanks, Matt > Best, > Hom Nath > > > > On Mon, Feb 8, 2016 at 2:27 PM, Barry Smith <bsm...@mcs.anl.gov> wrote: > > > > > > In this case you need to provide two pieces of information to the > PCFIELDSPLIT. What we call the "block size" or bs which is the number of > "basic fields" in the problem. For example if at each grid point you have > x-velocity, y-velocity, and pressure the block size is 3. The second is you > need to tell PCFIELDSPLIT what "basic fields" are in each split you want to > use. > > > > In your case you can do this with -pc_fieldsplit_block_size 4. > (Note if you use a DM to define your problem the PCFIELDSPLIT will > automatically get the bs from the DM so you do not need to provide it. > Similarly if you set a block size on your Mat the PCFIELDSPLIT will use > that so you don't need to set it). > > > > Then use -pc_fieldsplit_0_fields 0,1 to indicate your first split is > the 0 and 1 basic fields and -pc_fieldsplit_1_fields 2,3 to indicate the > second split is the 2 and 3 basic fields. > > (By default if you do not provide this then PCFIELDSPLIT will use bs > splits (4 in your case) one for each basic field. > > > > Barry > > > >> On Feb 8, 2016, at 11:20 AM, Hom Nath Gharti <hng.em...@gmail.com> > wrote: > >> > >> Hi Matt, Hi all, > >> > >> I am trying to have some feel for PCFIELDSPLIT testing on a very small > >> matrix (attached below). I have 4 degrees of freedoms. I use 4 > >> processors. I want to split 4 dofs into two fields each having two > >> dofs. I don't know whether this my be a problem for petsc. When I use > >> the command: > >> srun -n 4 ./xtestfs -ksp_monitor -ksp_view > >> > >> It runs fine. > >> > >> But when I want to use field split options using the command: > >> srun -n 4 ./xtestfs -ksp_monitor -ksp_type fgmres -pc_type fieldsplit > >> -fieldsplit_0_ksp_type gmres -fieldsplit_0_pc_type bjacobi > >> -fieldsplit_1_pc_type jacobi > >> > >> I get the following error message: > >> [0]PETSC ERROR: Petsc has generated inconsistent data > >> [0]PETSC ERROR: Unhandled case, must have at least two fields, not 1 > >> [0]PETSC ERROR: See > >> http://www.mcs.anl.gov/petsc/documentation/faq.html for trouble > >> shooting. > >> [0]PETSC ERROR: Petsc Release Version 3.6.3, Dec, 03, 2015 > >> [0]PETSC ERROR: /tigress/hgharti/test/testfs/./xtestfs on a > >> arch-linux2-c-debug named tiger-r11n1 by hgharti Mon Feb 8 11:40:03 > >> 2016 > >> [0]PETSC ERROR: Configure options > >> > --with-blas-lapack-dir=/opt/intel/composer_xe_2015.2.164/mkl/lib/intel64/ > >> --with-cc=mpicc --with-cxx=mpicxx --with-fc=mpif90 > >> --with-mpiexec=mpiexec --with-debugging=1 --download-scalapack > >> --download-mumps --download-pastix --download-superlu > >> --download-superlu_dist --download-metis --download-parmetis > >> --download-ptscotch --download-hypre > >> [0]PETSC ERROR: #1 PCFieldSplitSetDefaults() line 469 in > >> > /home/hgharti/lsoft/petsc-3.6.3-intel16-mpi5/src/ksp/pc/impls/fieldsplit/fieldsplit.c > >> [0]PETSC ERROR: #2 PCSetUp_FieldSplit() line 486 in > >> > /home/hgharti/lsoft/petsc-3.6.3-intel16-mpi5/src/ksp/pc/impls/fieldsplit/fieldsplit.c > >> [0]PETSC ERROR: #3 PCSetUp() line 983 in > >> > /home/hgharti/lsoft/petsc-3.6.3-intel16-mpi5/src/ksp/pc/interface/precon.c > >> [0]PETSC ERROR: #4 KSPSetUp() line 332 in > >> > /home/hgharti/lsoft/petsc-3.6.3-intel16-mpi5/src/ksp/ksp/interface/itfunc.c > >> [0]PETSC ERROR: #5 KSPSolve() line 546 in > >> > /home/hgharti/lsoft/petsc-3.6.3-intel16-mpi5/src/ksp/ksp/interface/itfunc.c > >> WARNING! There are options you set that were not used! > >> WARNING! could be spelling mistake, etc! > >> Option left: name:-fieldsplit_0_ksp_type value: gmres > >> Option left: name:-fieldsplit_0_pc_type value: bjacobi > >> Option left: name:-fieldsplit_1_pc_type value: jacobi > >> forrtl: error (76): Abort trap signal > >> > >> I tried several trials but could not fix the problem. Is it the > >> FORTRAN problem or am I doing something wrong? Any suggestions would > >> be greatly appreciated. > >> For your information I use: > >> PETSc-3.6.3 > >> intel-16.0 compiler > >> intel-mpi-5.1.1 > >> > >> Program > >> > >> > >> Best, > >> Hom Nath > >> > >> ! simple test case for PCFIELDSPLIT > >> ! > ----------------------------------------------------------------------- > >> program testfs > >> implicit none > >> #include "petsc/finclude/petscsys.h" > >> #include "petsc/finclude/petscvec.h" > >> #include "petsc/finclude/petscvec.h90" > >> #include "petsc/finclude/petscmat.h" > >> #include "petsc/finclude/petscksp.h" > >> #include "petsc/finclude/petscpc.h" > >> #include "petsc/finclude/petscviewer.h" > >> #include "petsc/finclude/petscviewer.h90" > >> > >> Vec x,b,u > >> Mat A > >> KSP ksp > >> PC pc > >> PetscErrorCode ierr > >> PetscBool flg > >> PetscInt errcode,its,maxitr,myrank,nproc > >> PetscInt nedof,nelmt,ndof,nenod,ngdof,fixdof,nsparse,neq > >> PetscInt,allocatable :: krow_sparse(:),ggdof_elmt(:,:), & > >> inode_elmt(:,:) > >> > >> PetscReal,allocatable :: storef(:,:),storekmat(:,:,:) > >> PetscInt gdof0(2),gdof1(2) > >> > >> ! initialize MPI > >> > >> call MPI_INIT(errcode) > >> if(errcode /= 0)stop 'ERROR: cannot initialize MPI!' > >> call MPI_COMM_RANK(MPI_COMM_WORLD,myrank,errcode) > >> if(errcode /= 0)stop 'ERROR: cannot find processor ID' > >> call MPI_COMM_SIZE(MPI_COMM_WORLD,nproc,errcode) > >> if(errcode /= 0)stop 'ERROR: cannot find number of processors!' > >> > >> ! define some parameters > >> ! 1D model with 4 elements. Each elements consits of 2 nodes. Node 0 is > fixed. > >> ! 1-----------2-----------3-----------4-----------5 > >> nelmt=1 ! per processor > >> nenod=2 ! number of nodes per element > >> nedof=2 ! number of DOFs per element > >> ndof=2 ! number of DOFs > >> ngdof=4 ! total number of global DOFs > >> fixdof=1 > >> > >> if(myrank==0)then > >> neq=1 ! number of equations > >> nsparse=1 > >> else > >> neq=2 ! number of equations > >> nsparse=4 > >> endif > >> > >> allocate(storef(nedof,nelmt),storekmat(nedof,nedof,nelmt), > & > >> ggdof_elmt(nedof,nelmt),inode_elmt(nenod,nelmt),krow_sparse(nsparse)) > >> > >> storef=0.0_8 ! local RHS vector > >> storekmat=0.0_8 ! local stiffness matrix > >> if(myrank==0)then > >> krow_sparse(:)=1 > >> storef(:,1)=(/12.5000_8, 12.5000_8/) > >> storekmat(:,:,1) = reshape((/1.2333333332_8, 0.0166666667_8, & > >> 0.0166666667_8, > 1.2333333332_8/),(/nedof,nedof/)) > >> inode_elmt(:,1) = (/1,2/) > >> ggdof_elmt(:,1) = (/0,1/) ! 0 because left node is fixed > >> elseif(myrank==1)then > >> krow_sparse(:)=(/1,1,2,2/) > >> storef(:,1)=(/12.5000_8, 12.5000_8/) > >> storekmat(:,:,1) = reshape((/1.2333333332_8, 0.0166666667_8, & > >> 0.0166666667_8, > 1.2333333332_8/),(/nedof,nedof/)) > >> inode_elmt(:,1) = (/1,2/) > >> ggdof_elmt(:,1) = (/1,2/) > >> elseif(myrank==2)then > >> krow_sparse(:)=(/1,1,2,2/) > >> storef(:,1)=(/12.5000_8, 12.5000_8/) > >> storekmat(:,:,1) = reshape((/1.2333333332_8, 0.0166666667_8, & > >> 0.0166666667_8, > 1.2333333332_8/),(/nedof,nedof/)) > >> inode_elmt(:,1) = (/1,2/) > >> ggdof_elmt(:,1) = (/2,3/) > >> elseif(myrank==3)then > >> krow_sparse(:)=(/1,1,2,2/) > >> storef(:,1)=(/12.5000_8, 22.5000_8/) > >> storekmat(:,:,1) = reshape((/1.2333333332_8, 0.0166666667_8, & > >> 0.0166666667_8, > 1.2333333332_8/),(/nedof,nedof/)) > >> inode_elmt(:,1) = (/1,2/) > >> ggdof_elmt(:,1) = (/3,4/) > >> endif > >> > >> call petsc_initialize > >> call petsc_set_matrix > >> call petsc_set_vector > >> call petsc_solve > >> call petsc_finalize > >> > >> ! deallocate > >> deallocate(storef,storekmat,ggdof_elmt,inode_elmt) > >> deallocate(krow_sparse) > >> contains > >> > !------------------------------------------------------------------------------- > >> subroutine petsc_initialize() > >> implicit none > >> PetscInt :: istart,iend > >> PetscInt :: nzeros_max,nzeros_min > >> PetscReal :: atol,dtol,rtol > >> PetscInt,allocatable :: nzeros(:) > >> IS gdof0_is,gdof1_is > >> > >> call PetscInitialize(PETSC_NULL_CHARACTER,ierr) > >> call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-n',ngdof,flg,ierr) > >> > >> ! count number of nonzeros per row > >> allocate(nzeros(neq)) > >> nzeros=0 > >> nzeros(krow_sparse)=nzeros(krow_sparse)+1 > >> nzeros_max=maxval(nzeros) > >> nzeros_min=minval(nzeros) > >> > >> ! create matrix > >> call MatCreate(PETSC_COMM_WORLD,A,ierr) > >> call MatSetType(A,MATMPIAIJ,ierr) > >> CHKERRQ(ierr) > >> call MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,ngdof,ngdof,ierr) > >> CHKERRQ(ierr) > >> ! preallocation > >> call > MatMPIAIJSetPreallocation(A,nzeros_max,PETSC_NULL_INTEGER,nzeros_max, & > >> PETSC_NULL_INTEGER,ierr) > >> CHKERRQ(ierr) > >> > >> call MatGetOwnershipRange(A,istart,iend,ierr) > >> CHKERRQ(ierr) > >> print*,'ownership:',myrank,istart,iend > >> ! create vector > >> call VecCreate(PETSC_COMM_WORLD,x,ierr) > >> call VecSetSizes(x,PETSC_DECIDE,ngdof,ierr) > >> call VecSetType(x,VECMPI,ierr) > >> call VecDuplicate(x,b,ierr) > >> call VecDuplicate(x,u,ierr) > >> > >> > !******************************************************************************* > >> ! get IS for split fields > >> gdof0=(/0,1/) > >> gdof1=(/2,3/) > >> call > ISCreateGeneral(PETSC_COMM_WORLD,2,gdof0,PETSC_COPY_VALUES,gdof0_is,ierr) > >> CHKERRQ(ierr) > >> call > ISCreateGeneral(PETSC_COMM_WORLD,2,gdof1,PETSC_COPY_VALUES,gdof1_is,ierr) > >> CHKERRQ(ierr) > >> > !******************************************************************************* > >> > >> ! Create linear solver context > >> call KSPCreate(PETSC_COMM_WORLD,ksp,ierr) > >> call KSPSetOperators(ksp,A,A,ierr) > >> > >> call KSPGetPC(ksp,pc,ierr) > >> CHKERRQ(ierr) > >> > >> > !******************************************************************************* > >> ! split PC > >> call PCFieldSplitSetIS(pc,"0",gdof0_is,ierr); > >> CHKERRQ(ierr) > >> call ISDestroy(gdof0_is,ierr) > >> CHKERRQ(ierr) > >> call PCFieldSplitSetIS(pc,"1",gdof1_is,ierr); > >> CHKERRQ(ierr) > >> call ISDestroy(gdof1_is,ierr) > >> CHKERRQ(ierr) > >> > !******************************************************************************* > >> > >> rtol = 1.0d-6 > >> atol = 1.0d-6 > >> dtol = 1.0d10 > >> maxitr = 1000 > >> call KSPSetTolerances(ksp,rtol,atol,dtol,maxitr,ierr) > >> CHKERRQ(ierr) > >> call KSPSetFromOptions(ksp,ierr) > >> end subroutine petsc_initialize > >> > !------------------------------------------------------------------------------- > >> > >> subroutine petsc_set_matrix() > >> > >> implicit none > >> integer :: i,i_elmt,j,ncount > >> integer :: ggdof(NEDOF),idof(NEDOF),igdof(NEDOF) > >> > >> ! Set and assemble matrix. > >> call MatZeroEntries(A,ierr) > >> CHKERRQ(ierr) > >> do i_elmt=1,nelmt > >> ggdof(:)=ggdof_elmt(:,i_elmt) > >> ggdof(:)=ggdof(:)-1 ! petsc index starts from 0 > >> ncount=0; idof=-1; igdof=-1 > >> do i=1,NEDOF > >> do j=1,NEDOF > >> if(ggdof(i).ge.0.and.ggdof(j).ge.0)then > >> call MatSetValues(A,1,ggdof(i),1,ggdof(j),storekmat(i,j,i_elmt), > >> & > >> ADD_VALUES,ierr) > >> CHKERRQ(ierr) > >> endif > >> enddo > >> enddo > >> enddo > >> call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr) > >> call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr) > >> if(myrank==0)print*,'matrix setting & assembly complete!' > >> > >> end subroutine petsc_set_matrix > >> > !------------------------------------------------------------------------------- > >> > >> subroutine petsc_set_vector() > >> implicit none > >> PetscScalar zero > >> integer :: i,i_elmt,ncount > >> integer :: ggdof(NEDOF),idof(NEDOF),igdof(NEDOF) > >> > >> ! set vector > >> zero=0.0_8 > >> call VecSet(b,zero,ierr) > >> do i_elmt=1,nelmt > >> ggdof(:)=ggdof_elmt(:,i_elmt) > >> ggdof(:)=ggdof(:)-1 ! petsc index starts from 0 > >> ncount=0; idof=-1; igdof=-1 > >> do i=1,NEDOF > >> if(ggdof(i).ge.0)then > >> call VecSetValues(b,1,ggdof(i),storef(i,i_elmt),ADD_VALUES,ierr); > >> CHKERRQ(ierr) > >> endif > >> enddo > >> enddo > >> > >> ! assemble vector > >> call VecAssemblyBegin(b,ierr) > >> call VecAssemblyEnd(b,ierr) > >> if(myrank==0)print*,'vector setting & assembly complete!' > >> > >> end subroutine petsc_set_vector > >> > !------------------------------------------------------------------------------- > >> > >> subroutine petsc_solve() > >> implicit none > >> PetscInt ireason > >> PetscScalar a_v(1) > >> PetscOffset a_i > >> PetscInt n > >> PetscReal,allocatable :: sdata(:) > >> > >> call VecGetSize(b,n,ierr) > >> CHKERRQ(ierr) > >> allocate(sdata(n)) > >> sdata=0.0_8 > >> call VecGetArray(b,a_v,a_i,ierr) > >> CHKERRQ(ierr) > >> sdata(1:n)=a_v(a_i+1:a_i+n) > >> call VecRestoreArray(b,a_v,a_i,ierr) > >> CHKERRQ(ierr) > >> print*,'rhs:',sdata > >> > >> call KSPSolve(ksp,b,x,ierr) > >> sdata=0.0_8 > >> call VecGetArray(x,a_v,a_i,ierr) > >> CHKERRQ(ierr) > >> sdata(1:n)=a_v(a_i+1:a_i+n) > >> call VecRestoreArray(b,a_v,a_i,ierr) > >> CHKERRQ(ierr) > >> print*,'solution:',sdata > >> > >> ! Check convergence > >> call KSPGetConvergedReason(ksp,ireason,ierr) > >> if(myrank==0)print*,'converges reason',myrank,ireason > >> call KSPGetIterationNumber(ksp,its,ierr) > >> if(myrank==0)print*,'Iterations:',its > >> deallocate(sdata) > >> end subroutine petsc_solve > >> > !------------------------------------------------------------------------------- > >> > >> subroutine petsc_finalize() > >> implicit none > >> > >> ! Free work space. > >> call VecDestroy(x,ierr) > >> call VecDestroy(u,ierr) > >> call VecDestroy(b,ierr) > >> call MatDestroy(A,ierr) > >> call KSPDestroy(ksp,ierr) > >> call PetscFinalize(ierr) > >> CHKERRQ(ierr) > >> > >> end subroutine petsc_finalize > >> > !------------------------------------------------------------------------------- > >> > >> end program testfs > >> > >> On Tue, Feb 2, 2016 at 4:54 PM, Matthew Knepley <knep...@gmail.com> > wrote: > >>> On Tue, Feb 2, 2016 at 3:11 PM, Hom Nath Gharti <hng.em...@gmail.com> > wrote: > >>>> > >>>> Thanks a lot. I will try. > >>> > >>> > >>> Also, if you send a small test case, we can run it and tell you the > problem. > >>> > >>> Matt > >>> > >>>> > >>>> Hom Nath > >>>> > >>>> On Tue, Feb 2, 2016 at 4:02 PM, Matthew Knepley <knep...@gmail.com> > wrote: > >>>>> On Tue, Feb 2, 2016 at 3:01 PM, Hom Nath Gharti <hng.em...@gmail.com > > > >>>>> wrote: > >>>>>> > >>>>>> Thanks again Matt. Unfortunately got the same errors with '0'. I > think > >>>>>> both are valid in Fortran. > >>>>> > >>>>> > >>>>> Then you will have to go in the debugger and see why its not > creating 4 > >>>>> splits, since this is exactly > >>>>> what it does in our tests. This is how all the SNES tests that I use > >>>>> work. I > >>>>> am sure its a Fortran > >>>>> problem. > >>>>> > >>>>> Thanks, > >>>>> > >>>>> Matt > >>>>> > >>>>>> > >>>>>> Hom Nath > >>>>>> > >>>>>> On Tue, Feb 2, 2016 at 3:42 PM, Matthew Knepley <knep...@gmail.com> > >>>>>> wrote: > >>>>>>> On Tue, Feb 2, 2016 at 2:35 PM, Hom Nath Gharti < > hng.em...@gmail.com> > >>>>>>> wrote: > >>>>>>>> > >>>>>>>> Thank you so much Matt. > >>>>>>>> > >>>>>>>> I now tried the following: > >>>>>>>> > >>>>>>>> ====================================================== > >>>>>>>> call KSPCreate(PETSC_COMM_WORLD,ksp,ierr) > >>>>>>>> call KSPGetPC(ksp,pc,ierr) > >>>>>>>> > >>>>>>>> call PCFieldSplitSetIS(pc,"0",gdofu_is,ierr); > >>>>>>>> call ISDestroy(gdofu_is,ierr) > >>>>>>>> call PCFieldSplitSetIS(pc,"1",gdofchi_is,ierr); > >>>>>>>> call ISDestroy(gdofchi_is,ierr) > >>>>>>>> call PCFieldSplitSetIS(pc,"2",gdofp_is,ierr); > >>>>>>>> call ISDestroy(gdofp_is,ierr) > >>>>>>>> call PCFieldSplitSetIS(pc,"3",gdofphi_is,ierr); > >>>>>>>> call ISDestroy(gdofphi_is,ierr) > >>>>>>>> > >>>>>>>> ! Amat changes in each time steps > >>>>>>>> call KSPSetOperators(ksp,Amat,Amat,ierr) !version >= 3.5.0 > >>>>>>>> > >>>>>>>> call KSPSolve(ksp,bvec,xvec,ierr) > >>>>>>>> ====================================================== > >>>>>>> > >>>>>>> > >>>>>>> I am guessing that "0" is not a valid string for your Fortran > >>>>>>> compiler. > >>>>>>> Are > >>>>>>> you sure > >>>>>>> its not single quotes '0'? > >>>>>>> > >>>>>>> Matt > >>>>>>> > >>>>>>>> > >>>>>>>> But I get the following error: > >>>>>>>> > >>>>>>>> [0]PETSC ERROR: --------------------- Error Message > >>>>>>>> -------------------------------------------------------------- > >>>>>>>> [0]PETSC ERROR: Petsc has generated inconsistent data > >>>>>>>> [0]PETSC ERROR: Unhandled case, must have at least two fields, > not 1 > >>>>>>>> [0]PETSC ERROR: See > >>>>>>>> http://www.mcs.anl.gov/petsc/documentation/faq.html for trouble > >>>>>>>> shooting. > >>>>>>>> [0]PETSC ERROR: Petsc Release Version 3.6.3, Dec, 03, 2015 > >>>>>>>> [0]PETSC ERROR: > >>>>>>>> /tigress/hgharti/gitwork/SPECFEM3D_GLOBEVSPP/./bin/xspecfem3D on a > >>>>>>>> arch-linux2-c-debug named tiger-r3c1n7 by hgharti Tue Feb 2 15: > >>>>>>>> 29:30 2016 > >>>>>>>> [0]PETSC ERROR: Configure options > >>>>>>>> > >>>>>>>> > >>>>>>>> > --with-blas-lapack-dir=/opt/intel/composer_xe_2015.2.164/mkl/lib/intel64/ > >>>>>>>> --with-cc=mpicc --with-cxx=mpicxx --wit > >>>>>>>> h-fc=mpif90 --with-mpiexec=mpiexec --with-debugging=1 > >>>>>>>> --download-scalapack --download-mumps --download-pastix > >>>>>>>> --download-superlu --download-superlu_dist --download-metis > >>>>>>>> --download-parmetis --download-ptscotch --download-hypre > >>>>>>>> [0]PETSC ERROR: #1 PCFieldSplitSetDefaults() line 469 in > >>>>>>>> > >>>>>>>> > >>>>>>>> > >>>>>>>> > /home/hgharti/lsoft/petsc-3.6.3-intel16-mpi5/src/ksp/pc/impls/fieldsplit/fieldsplit.c > >>>>>>>> [0]PETSC ERROR: #2 PCSetUp_FieldSplit() line 486 in > >>>>>>>> > >>>>>>>> > >>>>>>>> > >>>>>>>> > /home/hgharti/lsoft/petsc-3.6.3-intel16-mpi5/src/ksp/pc/impls/fieldsplit/fieldsplit.c > >>>>>>>> [0]PETSC ERROR: #3 PCSetUp() line 983 in > >>>>>>>> > >>>>>>>> > >>>>>>>> > /home/hgharti/lsoft/petsc-3.6.3-intel16-mpi5/src/ksp/pc/interface/precon.c > >>>>>>>> [0]PETSC ERROR: #4 KSPSetUp() line 332 in > >>>>>>>> > >>>>>>>> > >>>>>>>> > >>>>>>>> > /home/hgharti/lsoft/petsc-3.6.3-intel16-mpi5/src/ksp/ksp/interface/itfunc.c > >>>>>>>> [0]PETSC ERROR: #5 KSPSolve() line 546 in > >>>>>>>> > >>>>>>>> > >>>>>>>> > >>>>>>>> > /home/hgharti/lsoft/petsc-3.6.3-intel16-mpi5/src/ksp/ksp/interface/itfunc.c > >>>>>>>> forrtl: error (76): Abort trap signal > >>>>>>>> > >>>>>>>> Am I missing something? > >>>>>>>> > >>>>>>>> Thanks, > >>>>>>>> Hom Nath > >>>>>>>> > >>>>>>>> On Tue, Feb 2, 2016 at 3:24 PM, Matthew Knepley < > knep...@gmail.com> > >>>>>>>> wrote: > >>>>>>>>> On Tue, Feb 2, 2016 at 2:20 PM, Hom Nath Gharti > >>>>>>>>> <hng.em...@gmail.com> > >>>>>>>>> wrote: > >>>>>>>>>> > >>>>>>>>>> Hi Matt, hi all, > >>>>>>>>>> > >>>>>>>>>> I am trying to use PCFIELDSPLIT for my solver which consists of > 4 > >>>>>>>>>> different variables, namely, u (displacement vector), \chi > >>>>>>>>>> (displacement potential), p(pressure), and \phi (gravity > >>>>>>>>>> potential). > >>>>>>>>>> > >>>>>>>>>> My code segment looks like the following: > >>>>>>>>>> ====================================================== > >>>>>>>>>> call KSPCreate(PETSC_COMM_WORLD,ksp,ierr) > >>>>>>>>>> call KSPGetPC(ksp,pc,ierr) > >>>>>>>>>> nsplit=4 > >>>>>>>>>> call PCFieldSplitSetBlockSize(pc,nsplit,ierr); > >>>>>>>>> > >>>>>>>>> > >>>>>>>>> You do not need the statement above. > >>>>>>>>> > >>>>>>>>>> > >>>>>>>>>> call PCFieldSplitSetIS(pc,"0",gdofu_is,ierr); > >>>>>>>>>> call ISDestroy(gdofu_is,ierr) > >>>>>>>>>> call PCFieldSplitSetIS(pc,"1",gdofchi_is,ierr); > >>>>>>>>>> call ISDestroy(gdofchi_is,ierr) > >>>>>>>>>> call PCFieldSplitSetIS(pc,"2",gdofp_is,ierr); > >>>>>>>>>> call ISDestroy(gdofp_is,ierr) > >>>>>>>>>> call PCFieldSplitSetIS(pc,"3",gdofphi_is,ierr); > >>>>>>>>>> call ISDestroy(gdofphi_is,ierr) > >>>>>>>>>> > >>>>>>>>>> call PCFieldSplitGetSubKSP(pc,nsplit,subksp,ierr); > >>>>>>>>> > >>>>>>>>> > >>>>>>>>> These SetOperators() calls are wrong. I am not sure why you want > >>>>>>>>> them > >>>>>>>>> here. > >>>>>>>>> Also, that means you do not need the call above. > >>>>>>>>> > >>>>>>>>> Thanks, > >>>>>>>>> > >>>>>>>>> Matt > >>>>>>>>> > >>>>>>>>>> > >>>>>>>>>> call KSPSetOperators(subksp(1),Amat,Amat,ierr); > >>>>>>>>>> call KSPSetOperators(subksp(2),Amat,Amat,ierr); > >>>>>>>>>> call KSPSetOperators(subksp(3),Amat,Amat,ierr); > >>>>>>>>>> call KSPSetOperators(subksp(4),Amat,Amat,ierr); > >>>>>>>>>> > >>>>>>>>>> call KSPSolve(ksp,bvec,xvec,ierr) > >>>>>>>>>> ====================================================== > >>>>>>>>>> > >>>>>>>>>> But I am getting the following error: > >>>>>>>>>> [79]PETSC ERROR: Null argument, when expecting valid pointer > >>>>>>>>>> [79]PETSC ERROR: Null Object: Parameter # 1 > >>>>>>>>>> [79]PETSC ERROR: #1 KSPSetOperators() line 536 in > >>>>>>>>>> /home/hgharti/lsoft/petsc-3.6.3-intel16-mpi5/src/ksp/ksp/interf > >>>>>>>>>> > >>>>>>>>>> It looks like I am doing something wrong in "call > >>>>>>>>>> KSPSetOperators" > >>>>>>>>>> but > >>>>>>>>>> I could not figure that out. > >>>>>>>>>> > >>>>>>>>>> Could anybody help me to fix this problem? I looked into almost > >>>>>>>>>> all > >>>>>>>>>> related examples but I could not really figure out the correct > >>>>>>>>>> steps > >>>>>>>>>> after "call PCFieldSplitSetIS". > >>>>>>>>>> > >>>>>>>>>> Any help will be greatly appreciated. > >>>>>>>>>> > >>>>>>>>>> Best, > >>>>>>>>>> Hom nath > >>>>>>>>>> > >>>>>>>>>> On Sun, Jan 24, 2016 at 7:14 PM, Hom Nath Gharti > >>>>>>>>>> <hng.em...@gmail.com> > >>>>>>>>>> wrote: > >>>>>>>>>>> Thank you so much Matt! I will try. > >>>>>>>>>>> > >>>>>>>>>>> Hom Nath > >>>>>>>>>>> > >>>>>>>>>>> On Sun, Jan 24, 2016 at 6:26 AM, Matthew Knepley > >>>>>>>>>>> <knep...@gmail.com> > >>>>>>>>>>> wrote: > >>>>>>>>>>>> On Fri, Jan 22, 2016 at 2:19 PM, Hom Nath Gharti > >>>>>>>>>>>> <hng.em...@gmail.com> > >>>>>>>>>>>> wrote: > >>>>>>>>>>>>> > >>>>>>>>>>>>> Dear all, > >>>>>>>>>>>>> > >>>>>>>>>>>>> I am new to PcFieldSplit. > >>>>>>>>>>>>> > >>>>>>>>>>>>> I have a matrix formed using MATMPIAIJ. Is it possible to use > >>>>>>>>>>>>> PCFIELDSPLIT operations in this type of matrix? Or does it > >>>>>>>>>>>>> have > >>>>>>>>>>>>> to > >>>>>>>>>>>>> be > >>>>>>>>>>>>> MATMPIBIJ or MATNEST format? > >>>>>>>>>>>> > >>>>>>>>>>>> > >>>>>>>>>>>> Yes, you can split AIJ. > >>>>>>>>>>>> > >>>>>>>>>>>>> > >>>>>>>>>>>>> If possible for MATMPIAIJ, could anybody provide me a simple > >>>>>>>>>>>>> example > >>>>>>>>>>>>> or few steps? Variables in the equations are displacement > >>>>>>>>>>>>> vector, > >>>>>>>>>>>>> scalar potential and pressure. > >>>>>>>>>>>> > >>>>>>>>>>>> > >>>>>>>>>>>> If you do not have a collocated discretization, then you have > >>>>>>>>>>>> to > >>>>>>>>>>>> use > >>>>>>>>>>>> > >>>>>>>>>>>> > >>>>>>>>>>>> > >>>>>>>>>>>> > >>>>>>>>>>>> > >>>>>>>>>>>> > >>>>>>>>>>>> > http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/PC/PCFieldSplitSetIS.html > >>>>>>>>>>>> > >>>>>>>>>>>> Thanks, > >>>>>>>>>>>> > >>>>>>>>>>>> Matt > >>>>>>>>>>>> > >>>>>>>>>>>>> > >>>>>>>>>>>>> Thanks for help. > >>>>>>>>>>>>> > >>>>>>>>>>>>> Hom Nath > >>>>>>>>>>>> > >>>>>>>>>>>> > >>>>>>>>>>>> > >>>>>>>>>>>> > >>>>>>>>>>>> -- > >>>>>>>>>>>> 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 > >>>>>>>>> > >>>>>>>>> > >>>>>>>>> > >>>>>>>>> > >>>>>>>>> -- > >>>>>>>>> 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 > >>>>>>> > >>>>>>> > >>>>>>> > >>>>>>> > >>>>>>> -- > >>>>>>> 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 > >>>>> > >>>>> > >>>>> > >>>>> > >>>>> -- > >>>>> 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 > >>> > >>> > >>> > >>> > >>> -- > >>> 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 > > > -- 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