[petsc-users] Problem with fortran version of ex29 in ksp
TAY wee-beng
zonexo at gmail.com
Thu Apr 26 17:21:00 CDT 2012
Hi,
I managed to remove the out of range error by removing
*num = num + 1* from
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)
In C, num will give the number 5 since num start from 0.
In fortran, num starts from 1 so num = num + 1 before
MatSetValuesStencil must be removed.
I still can't get the same ans. I will try to figure it out. Just to
confirm, this subroutine:
call MatNullSpaceRemove(nullspace,b,PETSC_NULL,ierr)
is not added in fortran, although it is used in C.
Yours sincerely,
TAY wee-beng
On 26/4/2012 9:29 PM, Matthew Knepley wrote:
> On Thu, Apr 26, 2012 at 3:26 PM, TAY wee-beng <zonexo at gmail.com
> <mailto: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
>> <mailto: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 <mailto: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 <mailto: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/20120427/502a644b/attachment-0001.htm>
More information about the petsc-users
mailing list