[petsc-users] Problem with fortran version of ex29 in ksp
Matthew Knepley
knepley at gmail.com
Thu Apr 26 14:29:21 CDT 2012
On Thu, Apr 26, 2012 at 3:26 PM, TAY wee-beng <zonexo at gmail.com> wrote:
> Hi,
>
> Thanks for pointing out my mistake. I'm still learning about the code. I
> corrected but I still get the error:
>
I am sure its another simple bug. Go in with the debugger. When it says you
are inserting a new nonzero,
back up the stack to your routine and see what indices you are tying to
insert. If its not part of the stencil,
you have found the bug.
Matt
> [0]PETSC ERROR: --------------------- Error Message
> ----------------------------
> --------
> [0]PETSC ERROR: Argument out of range!
> [0]PETSC ERROR: New nonzero at (2,4) caused a malloc!
> [0]PETSC ERROR:
> ----------------------------------------------------------------
>
> boundary condition region:
>
> num = 1
>
> if (j/=0) then
>
> v(num) = -rho*HxdHy
>
> col(MatStencil_i,num) = i
>
> col(MatStencil_j,num) = j-1
>
> num = num + 1
>
> end if
>
> if (i/=0) then
>
> v(num) = -rho*HydHx
>
> col(MatStencil_i,num) = i-1
>
> col(MatStencil_j,num) = j
>
> num = num + 1
>
> end if
>
> if (i/=mx-1) then
>
> v(num) = -rho*HydHx
>
> col(MatStencil_i,num) = i+1
>
> col(MatStencil_j,num) = j
>
> num = num + 1
>
> end if
>
> if (j/=my-1) then
>
> v(num) = -rho*HxdHy
>
> col(MatStencil_i,num) = i
>
> col(MatStencil_j,num) = j+1
>
> num = num + 1
>
> end if
>
> v(num) = (num/2.0)*rho*(HxdHy + HydHx)
>
> col(MatStencil_i,num) = i
>
> col(MatStencil_j,num) = j+1
>
> num = num + 1
>
> call
> MatSetValuesStencil(jac,i1,row,num,col,v,INSERT_VALUES,ierr)
>
> Yours sincerely,
>
> TAY wee-beng
>
>
> On 26/4/2012 2:45 PM, Matthew Knepley wrote:
>
> 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
>
>
--
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/5e426494/attachment-0001.htm>
More information about the petsc-users
mailing list