[petsc-users] Problem with fortran version of ex29 in ksp

TAY wee-beng zonexo at gmail.com
Tue Apr 24 15:32:13 CDT 2012


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)

!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

             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

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20120424/c24ed305/attachment.htm>


More information about the petsc-users mailing list