[petsc-users] Problem with fortran version of ex29 in ksp
Matthew Knepley
knepley at gmail.com
Tue Apr 24 15:37:32 CDT 2012
On Tue, Apr 24, 2012 at 4:32 PM, TAY wee-beng <zonexo at gmail.com> wrote:
> Hi,
>
> I'm trying to rewrite ex29 (KSP^Laplacian, 2d) from c to fortran version.
> The fortran version is based on the template of ex22f, which is in 3d.
>
> I have attached the code below. I have 2 problems.
>
> 1. When i3 = -3, in dirichlet BC, I got the same answer as ex29 in c
> version.
>
> However, when I change i3 to -4, I got the following error:
>
> [0]PETSC ERROR: --------------------- Error Message
> ----------------------------
> --------
> [0]PETSC ERROR: Argument out of range!
> [0]PETSC ERROR: New nonzero at (9,1) caused a malloc!
> [0]PETSC ERROR:
> ----------------------------------------------------------------
> --------
> [0]PETSC ERROR: Petsc Development HG revision:
> fc934c1d84bc6ba8e2686702a8a99539d
> 50af122 HG Date: Mon Apr 23 11:39:32 2012 -0500
> [0]PETSC ERROR: See docs/changes/index.html for recent updates.
> [0]PETSC ERROR: See docs/faq.html for hints about trouble shooting.
> [0]PETSC ERROR: See docs/index.html for manual pages.
> [0]PETSC ERROR:
> ----------------------------------------------------------------
> --------
> [0]PETSC ERROR: c:\obj_tmp\ex29f\Debug\ex29f.exe on a petsc-3.2 named
> USER-PC by
> User Tue Apr 24 22:08:46 2012
> [0]PETSC ERROR: Libraries linked from
> /cygdrive/d/Lib/petsc-3.2-dev_win32_vs2008
> /lib
> [0]PETSC ERROR: Configure run at Mon Apr 23 21:45:20 2012
> [0]PETSC ERROR: Configure options --with-cc="win32fe cl"
> --with-fc="win32fe ifor
> t" --with-cxx="win32fe cl" --with-mpi-dir=/cygdrive/d/Lib/MPICH2/
> --download-f-b
> las-lapack=1 --prefix=/cygdrive/d/Lib/petsc-3.2-dev_win32_vs2008
> --with-debuggin
> g=1 --useThreads=0
> [0]PETSC ERROR:
> ----------------------------------------------------------------
> --------
> [0]PETSC ERROR: MatSetValues_SeqAIJ() line 331 in
> src/mat/impls/aij/seq/C:\temp\
> PETSC-~1\src\mat\impls\aij\seq\aij.c
> [0]PETSC ERROR: MatSetValues() line 1148 in
> src/mat/interface/C:\temp\PETSC-~1\s
> rc\mat\INTERF~1\matrix.c
> [0]PETSC ERROR: MatSetValuesLocal() line 2003 in
> src/mat/interface/C:\temp\PETSC
> -~1\src\mat\INTERF~1\matrix.c
> [0]PETSC ERROR: MatSetValuesStencil() line 1379 in
> src/mat/interface/C:\temp\PET
> SC-~1\src\mat\INTERF~1\matrix.c
> [0]PETSC ERROR: --------------------- Error Message
> ----------------------------
> --------
> [0]PETSC ERROR: Argument out of range!
> [0]PETSC ERROR: New nonzero at (10,2) caused a malloc!
>
> ...
>
> What did I do wrong?
>
>
> 2. When I wanted to use the neumann BC, following ex29, I added :
>
> In subroutine ComputeRHS
>
> call
> MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,nullspace,ierr)
>
> call MatNullSpaceRemove(nullspace,b,PETSC_NULL,ierr)
>
> call MatNullSpaceDestroy(nullspace,ierr)
>
> and in subroutine ComputeMatrix
>
> call
> MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,nullspace,ierr)
>
> call MatSetNullSpace(jac,nullspace,ierr)
>
> call MatNullSpaceDestroy(nullspace,ierr)
>
> When I compile, it says unresolved external symbol MatNullSpaceRemove.
> Removing MatNullSpaceRemove, I can compile but running gives the error:
>
> [0]PETSC ERROR: --------------------- Error Message
> ----------------------------
> --------
> [0]PETSC ERROR: Corrupt argument:
> see http://www.mcs.anl.gov/petsc/documentation/faq.html#valgrind!
> [0]PETSC ERROR: Invalid Pointer to Object: Parameter # 4!
> [0]PETSC ERROR:
> ----------------------------------------------------------------
> --------
> [0]PETSC ERROR: Petsc Development HG revision:
> fc934c1d84bc6ba8e2686702a8a99539d
> 50af122 HG Date: Mon Apr 23 11:39:32 2012 -0500
> [0]PETSC ERROR: See docs/changes/index.html for recent updates.
> [0]PETSC ERROR: See docs/faq.html for hints about trouble shooting.
> [0]PETSC ERROR: See docs/index.html for manual pages.
> [0]PETSC ERROR:
> ----------------------------------------------------------------
> --------
> [0]PETSC ERROR: c:\obj_tmp\ex29f\Debug\ex29f.exe on a petsc-3.2 named
> USER-PC by
> User Tue Apr 24 22:30:31 2012
> [0]PETSC ERROR: Libraries linked from
> /cygdrive/d/Lib/petsc-3.2-dev_win32_vs2008
> /lib
> [0]PETSC ERROR: Configure run at Mon Apr 23 21:45:20 2012
> [0]PETSC ERROR: Configure options --with-cc="win32fe cl"
> --with-fc="win32fe ifor
> t" --with-cxx="win32fe cl" --with-mpi-dir=/cygdrive/d/Lib/MPICH2/
> --download-f-b
> las-lapack=1 --prefix=/cygdrive/d/Lib/petsc-3.2-dev_win32_vs2008
> --with-debuggin
> g=1 --useThreads=0
> [0]PETSC ERROR:
> ----------------------------------------------------------------
> --------
> [0]PETSC ERROR: MatNullSpaceCreate() line 250 in
> src/mat/interface/C:\temp\PETSC
> -~1\src\mat\INTERF~1\matnull.c
> [0]PETSC ERROR: --------------------- Error Message
> ----------------------------
> --------
> [0]PETSC ERROR: Invalid argument!
> [0]PETSC ERROR: Wrong type of object: Parameter # 1!
> [0]PETSC ERROR:
> ----------------------------------------------------------------
> --------
> [0]PETSC ERROR: Petsc Development HG revision:
> fc934c1d84bc6ba8e2686702a8a99539d
> 50af122 HG Date: Mon Apr 23 11:39:32 2012 -0500
> [0]PETSC ERROR: See docs/changes/index.html for recent updates.
> [0]PETSC ERROR: See docs/faq.html for hints about trouble shooting.
> [0]PETSC ERROR: See docs/index.html for manual pages.
> [0]PETSC ERROR:
> ----------------------------------------------------------------
> --------
> [0]PETSC ERROR: c:\obj_tmp\ex29f\Debug\ex29f.exe on a petsc-3.2 named
> USER-PC by
> User Tue Apr 24 22:30:31 2012
> [0]PETSC ERROR: Libraries linked from
> /cygdrive/d/Lib/petsc-3.2-dev_win32_vs2008
> /lib
> [0]PETSC ERROR: Configure run at Mon Apr 23 21:45:20 2012
> [0]PETSC ERROR: Configure options --with-cc="win32fe cl"
> --with-fc="win32fe ifor
> t" --with-cxx="win32fe cl" --with-mpi-dir=/cygdrive/d/Lib/MPICH2/
> --download-f-b
> las-lapack=1 --prefix=/cygdrive/d/Lib/petsc-3.2-dev_win32_vs2008
> --with-debuggin
> g=1 --useThreads=0
> [0]PETSC ERROR:
> ----------------------------------------------------------------
> --------
> [0]PETSC ERROR: MatNullSpaceDestroy() line 305 in
> src/mat/interface/C:\temp\PETS
> C-~1\src\mat\INTERF~1\matnull.c
> [0]PETSC ERROR: ourdmfunction() line 30 in
> src/dm/interface/ftn-custom/C:\temp\P
> ETSC-~1\src\dm\INTERF~1\FTN-CU~1\zdmf.c
> [0]PETSC ERROR: DMComputeFunction() line 1783 in
> src/dm/interface/C:\temp\PETSC-
> ~1\src\dm\INTERF~1\dm.c
> [0]PETSC ERROR: KSPSetUp() line 218 in
> src/ksp/ksp/interface/C:\temp\PETSC-~1\sr
> c\ksp\ksp\INTERF~1\itfunc.c
> [0]PETSC ERROR: --------------------- Error Message
> ----------------------------
> --------
> [0]PETSC ERROR: Corrupt argument:
> see http://www.mcs.anl.gov/petsc/documentation/faq.html#valgrind!
> [0]PETSC ERROR: Invalid Pointer to Object: Parameter # 4!
> [0]PETSC ERROR:
> ----------------------------------------------------------------
> --------
> [0]PETSC ERROR: Petsc Development HG revision:
> fc934c1d84bc6ba8e2686702a8a99539d
> 50af122 HG Date: Mon Apr 23 11:39:32 2012 -0500
> [0]PETSC ERROR: See docs/changes/index.html for recent updates.
> [0]PETSC ERROR: See docs/faq.html for hints about trouble shooting.
> [0]PETSC ERROR: See docs/index.html for manual pages.
> [0]PETSC ERROR:
> ----------------------------------------------------------------
> --------
> [0]PETSC ERROR: c:\obj_tmp\ex29f\Debug\ex29f.exe on a petsc-3.2 named
> USER-PC by
> User Tue Apr 24 22:30:31 2012
> [0]PETSC ERROR: Libraries linked from
> /cygdrive/d/Lib/petsc-3.2-dev_win32_vs2008
> /lib
> [0]PETSC ERROR: Configure run at Mon Apr 23 21:45:20 2012
> [0]PETSC ERROR: Configure options --with-cc="win32fe cl"
> --with-fc="win32fe ifor
> t" --with-cxx="win32fe cl" --with-mpi-dir=/cygdrive/d/Lib/MPICH2/
> --download-f-b
> las-lapack=1 --prefix=/cygdrive/d/Lib/petsc-3.2-dev_win32_vs2008
> --with-debuggin
> g=1 --useThreads=0
> [0]PETSC ERROR:
> ----------------------------------------------------------------
> --------
> [0]PETSC ERROR: MatNullSpaceCreate() line 250 in
> src/mat/interface/C:\temp\PETSC
> -~1\src\mat\INTERF~1\matnull.c
> [0]PETSC ERROR: --------------------- Error Message
> ----------------------------
> --------
> [0]PETSC ERROR: Invalid argument!
> [0]PETSC ERROR: Wrong type of object: Parameter # 1!
> [0]PETSC ERROR:
> ----------------------------------------------------------------
> --------
> [0]PETSC ERROR: Petsc Development HG revision:
> fc934c1d84bc6ba8e2686702a8a99539d
> 50af122 HG Date: Mon Apr 23 11:39:32 2012 -0500
> [0]PETSC ERROR: See docs/changes/index.html for recent updates.
> [0]PETSC ERROR: See docs/faq.html for hints about trouble shooting.
> [0]PETSC ERROR: See docs/index.html for manual pages.
> [0]PETSC ERROR:
> ----------------------------------------------------------------
> --------
> [0]PETSC ERROR: c:\obj_tmp\ex29f\Debug\ex29f.exe on a petsc-3.2 named
> USER-PC by
> User Tue Apr 24 22:30:31 2012
> [0]PETSC ERROR: Libraries linked from
> /cygdrive/d/Lib/petsc-3.2-dev_win32_vs2008
> /lib
> [0]PETSC ERROR: Configure run at Mon Apr 23 21:45:20 2012
> [0]PETSC ERROR: Configure options --with-cc="win32fe cl"
> --with-fc="win32fe ifor
> t" --with-cxx="win32fe cl" --with-mpi-dir=/cygdrive/d/Lib/MPICH2/
> --download-f-b
> las-lapack=1 --prefix=/cygdrive/d/Lib/petsc-3.2-dev_win32_vs2008
> --with-debuggin
> g=1 --useThreads=0
> [0]PETSC ERROR:
> ----------------------------------------------------------------
> --------
> [0]PETSC ERROR: MatNullSpaceDestroy() line 305 in
> src/mat/interface/C:\temp\PETS
> C-~1\src\mat\INTERF~1\matnull.c
> [0]PETSC ERROR: ourdmfunction() line 30 in
> src/dm/interface/ftn-custom/C:\temp\P
> ETSC-~1\src\dm\INTERF~1\FTN-CU~1\zdmf.c
> [0]PETSC ERROR: DMComputeFunction() line 1783 in
> src/dm/interface/C:\temp\PETSC-
> ~1\src\dm\INTERF~1\dm.c
> [0]PETSC ERROR: KSPSetUp() line 218 in
> src/ksp/ksp/interface/C:\temp\PETSC-~1\sr
> c\ksp\ksp\INTERF~1\itfunc.c
> [0]PETSC ERROR: KSPSolve() line 402 in
> src/ksp/ksp/interface/C:\temp\PETSC-~1\sr
> c\ksp\ksp\INTERF~1\itfunc.c
> Vector Object:Vec_0000000084000000_0 1 MPI processes
> type: mpi
> Process [0]
>
> What's wrong?
>
>
> *program ex29f*
>
> implicit none
>
> ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> ! Include files
> ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> !
> ! petscsys.h - base PETSc routines petscvec.h - vectors
> ! petscmat.h - matrices
> ! petscksp.h - Krylov subspace methods petscpc.h - preconditioners
>
> #include "finclude/petsc.h90"
>
> PetscErrorCode ierr
>
> DM da
>
> KSP ksp
>
> Vec x
>
> external ComputeRHS,ComputeMatrix
>
> PetscInt i1,i3
>
> call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
>
> i3 = -4
>
> i1 = 1
>
> call KSPCreate(MPI_COMM_WORLD,ksp,ierr)
>
> call
> DMDACreate2d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,i3,i3,PETSC_DECIDE,PETSC_DECIDE,i1,i1,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,da,ierr)
> call DMSetFunction(da,ComputeRHS,ierr)
> call DMSetJacobian(da,ComputeMatrix,ierr)
> call KSPSetDM(ksp,da,ierr)
>
> call KSPSetFromOptions(ksp,ierr)
> call KSPSetUp(ksp,ierr)
> call KSPSolve(ksp,PETSC_NULL_OBJECT,PETSC_NULL_OBJECT,ierr)
> call KSPGetSolution(ksp,x,ierr)
> call VecView(x,PETSC_VIEWER_STDOUT_WORLD,ierr)
> call KSPDestroy(ksp,ierr)
> call DMDestroy(da,ierr)
> call PetscFinalize(ierr)
>
> end program ex29f
>
> subroutine ComputeRHS(da,x,b,ierr)
>
> implicit none
>
> #include "finclude/petsc.h90"
>
> PetscErrorCode ierr
> PetscInt i,j,mx,my,xm,ym,xs,ys
> PetscScalar h,nu,rho
> PetscScalar Hx,Hy
> PetscScalar,pointer :: array(:,:)
> Vec x,b
> DM da
> MatNullSpace nullspace !>neumann BC
>
> nu = 0.1
>
> call
> DMDAGetInfo(da,PETSC_NULL_INTEGER,mx,my,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr)
>
> Hx = 1.d0 / (mx-1)
>
> Hy = 1.d0 / (my-1)
>
> call
> DMDAGetCorners(da,xs,ys,PETSC_NULL_INTEGER,xm,ym,PETSC_NULL_INTEGER,ierr)
>
> call DMDAVecGetArrayF90(da,b,array,ierr)
>
> do j = ys,ys+ym-1
>
> do i = xs,xs+xm-1
>
> array(i,j) = exp(-(i*Hx)*(i*Hx)/nu)*exp(-(j*Hy)*(j*Hy)/nu)*Hx*Hy
>
> end do
>
> end do
>
> call DMDAVecRestoreArrayF90(da,b,array,ierr)
>
> call VecAssemblyBegin(b,ierr)
>
> call VecAssemblyEnd(b,ierr)
>
> !call
> MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,nullspace,ierr)
>
4th argument should be PETSC_NULL_OBJECT.
Matt
> !call MatNullSpaceRemove(nullspace,b,PETSC_NULL,ierr)
>
> !call MatNullSpaceDestroy(nullspace,ierr)
>
> end subroutine ComputeRHS
>
>
> subroutine ComputeMatrix(da,x,JJ,jac,str,ierr)
>
> implicit none
>
> #include "finclude/petsc.h90"
>
> Mat jac,JJ
>
> PetscErrorCode ierr
>
> DM da
>
> PetscInt i,j,k,mx,my,xm
>
> PetscInt ym,xs,ys,i1,i5
>
> PetscScalar v(5),Hx,Hy,rho,centerRho
>
> PetscScalar HydHx
>
> PetscScalar HxdHy
>
> MatStencil row(4),col(4,5)
>
> Vec x
>
> MatStructure str
>
> MatNullSpace nullspace !>neumann BC
>
> rho = 1.0
>
> i1 = 1
>
> i5 = 5
>
> centerRho = rho
>
> call
> DMDAGetInfo(da,PETSC_NULL_INTEGER,mx,my,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr)
>
> Hx = 1.d0 / (mx-1)
>
> Hy = 1.d0 / (my-1)
>
> HxdHy = Hx/Hy
>
> HydHx = Hy/Hx
>
> call
> DMDAGetCorners(da,xs,ys,PETSC_NULL_INTEGER,xm,ym,PETSC_NULL_INTEGER,ierr)
>
>
> do j=ys,ys+ym-1
>
> do i=xs,xs+xm-1
>
> row(MatStencil_i) = i
>
> row(MatStencil_j) = j
>
> call ComputeRho(i,j,mx,my,centerRho,rho)
>
> if (i.eq.0 .or. j.eq.0 .or. i.eq.mx-1 .or. j.eq.my-1) then
>
> v(1) = 2.0*rho*(HxdHy + HydHx)
>
> call
> MatSetValuesStencil(jac,i1,row,i1,row,v,INSERT_VALUES,ierr)
>
> else
>
> v(1) = -rho*HxdHy
>
> col(MatStencil_i,1) = i
>
> col(MatStencil_j,2) = j-1
>
Cut and paste error above.
Matt
> v(2) = -rho*HydHx
>
> col(MatStencil_i,2) = i-1
>
> col(MatStencil_j,2) = j
>
> v(3) = 2.0*rho*(HxdHy + HydHx)
>
> col(MatStencil_i,3) = i
>
> col(MatStencil_j,3) = j
>
> v(4) = -rho*HydHx
>
> col(MatStencil_i,4) = i+1
>
> col(MatStencil_j,4) = j
>
> v(5) = -rho*HxdHy
>
> col(MatStencil_i,5) = i
>
> col(MatStencil_j,5) = j+1
>
> call
> MatSetValuesStencil(jac,i1,row,i5,col,v,INSERT_VALUES,ierr)
>
> end if
>
> end do
>
> end do
>
> call MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr)
>
> call MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr)
>
> !call
> MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,nullspace,ierr)
>
> !call MatSetNullSpace(jac,nullspace,ierr)
>
> !call MatNullSpaceDestroy(nullspace,ierr)
>
> end subroutine ComputeMatrix
>
> subroutine ComputeRho(i,j,mx,my,centerRho,rho)
>
> PetscInt i,j,mx,my
>
> PetscScalar rho,centerRho
>
> if ((i > mx/3.0) .and. (i < 2.0*mx/3.0) .and. (j > my/3.0) .and. (j <
> 2.0*my/3.0)) then
>
> rho = centerRho
>
> else
>
> rho = 1.0
>
> end if
>
>
> end subroutine ComputeRho
>
> --
> Yours sincerely,
>
> TAY wee-beng
>
>
--
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/20120424/4795406c/attachment-0001.htm>
More information about the petsc-users
mailing list