[petsc-users] Problem with fortran version of ex29 in ksp
TAY wee-beng
zonexo at gmail.com
Tue May 1 16:48:19 CDT 2012
Hi,
Do you mean my method is wrong?
I am following the template of ex22f,
where the variables are declared as :
PetscScalar v(5)
MatStencil row(4),col(4,5)
Hence,
for the neumann BC
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-1)/2.0)*rho*(HxdHy + HydHx)
print *, v
col(MatStencil_i,num) = i
col(MatStencil_j,num) = j
!num = num + 1
call
MatSetValuesStencil(jac,i1,row,num,col,v,INSERT_VALUES,ierr)
I do not get any more out of range error. However,my ans is still
different from that of ex29 in C.
Yours sincerely,
TAY wee-beng
On 27/4/2012 1:22 AM, Barry Smith wrote:
> From the manual page for MatSetValuesStencil():
>
> MatSetValuesStencil() uses 0-based row and column numbers in Fortran
> as well as in C.
>
>
> On Apr 26, 2012, at 5:21 PM, TAY wee-beng wrote:
>
>> 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> 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:
>>>>>> seehttp://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:
>>>>>> seehttp://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
More information about the petsc-users
mailing list