[petsc-users] Problem with fortran version of ex29 in ksp
Matthew Knepley
knepley at gmail.com
Thu Apr 26 07:45:22 CDT 2012
On Thu, Apr 26, 2012 at 8:27 AM, TAY wee-beng <zonexo at gmail.com> wrote:
> Hi,
>
> Is there an answer to the question below?
> I still can't get it working.
>
> Thanks.
>
>
> Hi,
>
> There's no error now but I didn't get the same ans as ex29. I believe the
> error lies in the boundary evaluation of v.
>
> In c, the code is :
>
> else if (user->bcType == NEUMANN) {
> num = 0;
> if (j!=0) {
> v[num] = -rho*HxdHy; col[num].i = i; col[num].j
> = j-1;
>
^^^ See
this?
> num++;
> }
> if (i!=0) {
> v[num] = -rho*HydHx; col[num].i = i-1; col[num].j
> = j;
> num++;
> }
> if (i!=mx-1) {
> v[num] = -rho*HydHx; col[num].i = i+1; col[num].j
> = j;
> num++;
> }
> if (j!=my-1) {
> v[num] = -rho*HxdHy; col[num].i = i; col[num].j
> = j+1;
> num++;
> }
> v[num] = (num/2.0)*rho*(HxdHy + HydHx); col[num].i = i;
> col[num].j = j;
> num++;
> ierr =
> MatSetValuesStencil(jac,1,&row,num,col,v,INSERT_VALUES);CHKERRQ(ierr);
>
> In fortran, I have changed to :
>
> else if (bc_cond == 2) then
>
> num = 1
>
> if (j/=0) then
>
> v(num) = -rho*HxdHy
>
> col(MatStencil_i,1) = i
>
^^ Why is this 1 in stead of
'num'?
Matt
col(MatStencil_j,2) = j-1
>
> num = num + 1
>
> end if
>
> if (i/=0) then
>
> v(num) = -rho*HydHx
>
> col(MatStencil_i,1) = i-1
>
> col(MatStencil_j,2) = j
>
> num = num + 1
>
> end if
>
> if (i/=mx-1) then
>
> v(num) = -rho*HydHx
>
> col(MatStencil_i,1) = i+1
>
> col(MatStencil_j,2) = j
>
> num = num + 1
>
> end if
>
> if (j/=my-1) then
>
> v(num) = -rho*HxdHy
>
> col(MatStencil_i,1) = i
>
> col(MatStencil_j,2) = j+1
>
> num = num + 1
>
> end if
>
> v(num) = (num/2.0)*rho*(HxdHy + HydHx)
>
> col(MatStencil_i,1) = i
>
> col(MatStencil_j,2) = j+1
>
> num = num + 1
>
> call
> MatSetValuesStencil(jac,i1,row,num,col,v,INSERT_VALUES,ierr)
>
> However, I got the error:
>
> [0]PETSC ERROR: --------------------- Error Message
> ---------------------------
> --------
> [0]PETSC ERROR: Argument out of range!
> [0]PETSC ERROR: New nonzero at (1,3) caused a malloc!
> [0]PETSC ERROR:
> ---------------------------------------------------------------
>
> How such is be defined in Fortran?
>
> Thank you
>
> Yours sincerely,
>
> TAY wee-beng
>
>
> On 25/4/2012 12:06 AM, Matthew Knepley wrote:
>
> On Tue, Apr 24, 2012 at 4:57 PM, TAY wee-beng <zonexo at gmail.com> wrote:
>
>> Hi,
>>
>> I still got the same error the moment this subroutine is called, after
>> changing to :
>>
>> call
>> MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,PETSC_NULL_INTEGER,PETSC_NULL_OBJECT,nullspace,ierr)
>>
>
> Argument 3 should be 0.
>
> Matt
>
>
>> [0]PETSC ERROR: Null argument, when expecting valid pointer!
>> [0]PETSC ERROR: Null 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:54:50 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
>>
>> Yours sincerely,
>>
>> TAY wee-beng
>>
>>
>> On 24/4/2012 10:37 PM, Matthew Knepley wrote:
>>
>> 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
>>
>>
>
>
> --
> 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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20120426/4a2ae495/attachment-0001.htm>
More information about the petsc-users
mailing list