[petsc-users] Fwd: Re: Problem with fortran version of ex29 in ksp
TAY wee-beng
zonexo at gmail.com
Sun Jun 3 11:13:15 CDT 2012
Hi,
I have update the ex29f.F90 since there's some problems with MPI in the
old version. Hopefully the DMDAVecGetArrayF90 and DMDAVecRestoreArrayF90
can work in linux.
Yours sincerely,
TAY wee-beng
On 11/5/2012 4:04 PM, John Mousel wrote:
> Include both header files.
>
> #include <finclude/petscdmda.h>
> #include <finclude/petscdmda.h90>
>
> The h90 header only includes the interface definitions for Fortran.
>
> John
>
> On Fri, May 11, 2012 at 8:51 AM, TAY wee-beng <zonexo at gmail.com
> <mailto:zonexo at gmail.com>> wrote:
>
>
>
> On 11/5/2012 3:30 PM, Matthew Knepley wrote:
>> On Fri, May 11, 2012 at 9:24 AM, TAY wee-beng <zonexo at gmail.com
>> <mailto:zonexo at gmail.com>> wrote:
>>
>> Hi,
>>
>> I have been using the GUI environment to do debugging so I am
>> a bit reluctant to learn Valgrind as its outputs seems a bit
>> daunting. But I guess John is right. I've been spending
>> these few days learning bit by bit.
>>
>> I realised that the error occurs in computerhs, at:
>>
>>
>> I bet this is a beautiful Fortranism. Do you include the F90
>> header file with the interface definition?
>> If not, Fortran just craps out like this. I can't stress enough
>> how much time would be saved by
>> switching languages to something with at least a modicum of error
>> checking.
>
> I initially used:
>
> #include "finclude/petsc.h90"
>
> Compilation and linking was fine in Linux and vs2008.
>
> Now I changed it to what ex22.F was using :
>
> #include <finclude/petscsys.h>
> #include <finclude/petscvec.h>
> #include <finclude/petscmat.h>
> #include <finclude/petscpc.h>
> #include <finclude/petscksp.h>
> #include <finclude/petscdmda.h>
>
> Compiling was ok but linking failed in Linux and VS2008:
>
> undefined reference to `dmdavecgetarrayf90_'
>
> I tried changing #include <finclude/petscdmda.h> to #include
> <finclude/petscdmda.h90> and everything was ok in VS2008 again,
> giving the right answers.
>
> However, in Linux, I got the following error:
>
> [wtay at hpc12:tutorials]$ /opt/openmpi-1.5.3/bin/mpif90 -c -fPIC
> -g -I/home/wtay/Lib/petsc-3.2-dev_shared_debug/include
> -I/home/wtay/Lib/petsc-3.2-dev_shared_debug/include
> -I/opt/openmpi-1.5.3/include -o ex29f.o ex29f.F90
> /home/wtay/Lib/petsc-3.2-dev_shared_debug/include/finclude/ftn-custom/petscdmda.h90(25):
> error #5082: Syntax error, found '::' when expecting one of: ( % :
> . = =>
> DMDABoundaryType :: pt
> -------------------------^
> /home/wtay/Lib/petsc-3.2-dev_shared_debug/include/finclude/ftn-custom/petscdmda.h90(26):
> error #5082: Syntax error, found '::' when expecting one of: ( % :
> . = =>
> DMDAStencilType :: st
> -------------------------^
> /home/wtay/Lib/petsc-3.2-dev_shared_debug/include/finclude/ftn-custom/petscdmda.h90(25):
> error #6590: This statement is not permitted as a statement within
> a derived-type-def
> DMDABoundaryType :: pt
> --------^
> /home/wtay/Lib/petsc-3.2-dev_shared_debug/include/finclude/ftn-custom/petscdmda.h90(26):
> error #6590: This statement is not permitted as a statement within
> a derived-type-def
> DMDAStencilType :: st
> --------^
> ex29f.F90(68): error #6404: This name does not have a type, and
> must have an explicit type. [DMDA_BOUNDARY_NONE]
> 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)
> -----------------------------------^
> ex29f.F90(68): error #6404: This name does not have a type, and
> must have an explicit type. [DMDA_STENCIL_STAR]
> 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)
> -------------------------------------------------------------------------^
> /home/wtay/Lib/petsc-3.2-dev_shared_debug/include/finclude/ftn-custom/petscdmda.h90(25):
> error #5082: Syntax error, found '::' when expecting one of: ( % :
> . = =>
> DMDABoundaryType :: pt
> -------------------------^
> /home/wtay/Lib/petsc-3.2-dev_shared_debug/include/finclude/ftn-custom/petscdmda.h90(26):
> error #5082: Syntax error, found '::' when expecting one of: ( % :
> . = =>
> DMDAStencilType :: st
> -------------------------^
> /home/wtay/Lib/petsc-3.2-dev_shared_debug/include/finclude/ftn-custom/petscdmda.h90(25):
> error #6590: This statement is not permitted as a statement within
> a derived-type-def
> DMDABoundaryType :: pt
> --------^
> /home/wtay/Lib/petsc-3.2-dev_shared_debug/include/finclude/ftn-custom/petscdmda.h90(26):
> error #6590: This statement is not permitted as a statement within
> a derived-type-def
> DMDAStencilType :: st
> --------^
> /home/wtay/Lib/petsc-3.2-dev_shared_debug/include/finclude/ftn-custom/petscdmda.h90(25):
> error #5082: Syntax error, found '::' when expecting one of: ( % :
> . = =>
> DMDABoundaryType :: pt
> -------------------------^
> /home/wtay/Lib/petsc-3.2-dev_shared_debug/include/finclude/ftn-custom/petscdmda.h90(26):
> error #5082: Syntax error, found '::' when expecting one of: ( % :
> . = =>
> DMDAStencilType :: st
> -------------------------^
> /home/wtay/Lib/petsc-3.2-dev_shared_debug/include/finclude/ftn-custom/petscdmda.h90(25):
> error #6590: This statement is not permitted as a statement within
> a derived-type-def
> DMDABoundaryType :: pt
> --------^
> /home/wtay/Lib/petsc-3.2-dev_shared_debug/include/finclude/ftn-custom/petscdmda.h90(26):
> error #6590: This statement is not permitted as a statement within
> a derived-type-def
> DMDAStencilType :: st
>
> Is there some errors in petscdmda.h90?
>
>
>
>>
>> Matt
>>
>> *call DMDAVecRestoreArrayF90(da,b,array,ierr)*
>>
>> ==27464== Invalid write of size 8
>> ==27464== at 0x402835: computerhs_ (ex29f.F90:119)
>> ==27464== Address 0xfffffffffffffef0 is not stack'd,
>> malloc'd or (recently) free'd
>> ==27464==
>> [0]PETSC ERROR:
>> ------------------------------------------------------------------------
>> [0]PETSC ERROR: Caught signal number 11 SEGV: Segmentation
>> Violation, probably memory access out of range
>> [0]PETSC ERROR: Try option -start_in_debugger or
>> -on_error_attach_debugger
>> [0]PETSC ERROR: or see
>> http://www.mcs.anl.gov/petsc/documentation/faq.html#valgrind[0]PETSC
>> ERROR: or try http://valgrind.org on GNU/linux and Apple Mac
>> OS X to find memory corruption errors
>> [0]PETSC ERROR: likely location of problem given in stack below
>> [0]PETSC ERROR: --------------------- Stack Frames
>> ------------------------------------
>> [0]PETSC ERROR: Note: The EXACT line numbers in the stack are
>> not available,
>> [0]PETSC ERROR: INSTEAD the line number of the start of
>> the function
>> [0]PETSC ERROR: is given.
>> [0]PETSC ERROR: [0] DM user function line 0 unknownunknown
>> [0]PETSC ERROR: [0] DMComputeFunction line 2085
>> /home/wtay/Codes/petsc-dev/src/dm/interface/dm.c
>> [0]PETSC ERROR: [0] KSPSetUp line 182
>> /home/wtay/Codes/petsc-dev/src/ksp/ksp/interface/itfunc.c
>> [0]PETSC ERROR: --------------------- Error Message
>> ------------------------------------
>> [0]PETSC ERROR: Signal received!
>>
>> I have checked that "array" 's values are correct. This
>> statement executed without problems in VS2008. If I replace
>> the above with something like:
>>
>> *call VecSet(b,Hy,ierr)*
>>
>> Everything is fine in Linux.
>>
>> Is there something wrong with *DMDAVecRestoreArrayF90*?
>>
>> My code in the area is:
>>
>> 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)
>>
>>
>> Yours sincerely,
>>
>> TAY wee-beng
>>
>>
>> On 8/5/2012 9:41 PM, John Mousel wrote:
>>> TAY wee-bing,
>>>
>>> If you want to be a programmer that writes interesting and
>>> reliable code, you need to be willing to use the tools of
>>> the trade. I can't think of a bigger time-saver than
>>> Valgrind. I would suggest that you learn to use it and use
>>> it a lot. I bet it will lead you to the root of your problem
>>> pretty quickly.
>>>
>>> John
>>>
>>> On Tue, May 8, 2012 at 2:17 PM, TAY wee-beng
>>> <zonexo at gmail.com <mailto:zonexo at gmail.com>> wrote:
>>>
>>> Hi,
>>>
>>> I compiled and run my code under visual studio 2008 with
>>> intel fortran. Everything works ok.
>>>
>>> However, when I tried to run the code in linux, I got
>>> the error as below. The error happens when
>>> KSPSetUp(ksp,ierr) is called.
>>>
>>> However, I am not able to print VecView or MatView to
>>> view if there's any errors. Is there any recommendation
>>> for debugging? I hope I do not need to valgrind if possible.
>>>
>>> [0]PETSC ERROR:
>>> ------------------------------------------------------------------------
>>> [0]PETSC ERROR: Caught signal number 11 SEGV:
>>> Segmentation Violation, probably memory access out of range
>>> [0]PETSC ERROR: Try option -start_in_debugger or
>>> -on_error_attach_debugger
>>> [0]PETSC ERROR: or see
>>> http://www.mcs.anl.gov/petsc/documentation/faq.html#valgrind[0]PETSC
>>> ERROR: or try http://valgrind.org on GNU/linux and Apple
>>> Mac OS X to find memory corruption errors
>>> [0]PETSC ERROR: likely location of problem given in
>>> stack below
>>> [0]PETSC ERROR: --------------------- Stack Frames
>>> ------------------------------------
>>> [0]PETSC ERROR: Note: The EXACT line numbers in the
>>> stack are not available,
>>> [0]PETSC ERROR: INSTEAD the line number of the
>>> start of the function
>>> [0]PETSC ERROR: is given.
>>> [0]PETSC ERROR: [0] DM user function line 0 unknownunknown
>>> [0]PETSC ERROR: [0] DMComputeFunction line 2085
>>> /home/wtay/Codes/petsc-dev/src/dm/interface/dm.c
>>> [0]PETSC ERROR: [0] KSPSetUp line 182
>>> /home/wtay/Codes/petsc-dev/src/ksp/ksp/interface/itfunc.c
>>> [0]PETSC ERROR: --------------------- Error Message
>>> ------------------------------------
>>> [0]PETSC ERROR: Signal received!
>>> [0]PETSC ERROR:
>>> ------------------------------------------------------------------------
>>> [0]PETSC ERROR: Petsc Development HG revision:
>>> 7ecdd63ec420b1659b960e65d96e822c5ac1a968 HG Date: Mon
>>> May 07 21:42:26 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: ./ex29f on a petsc-3.2 named hpc12 by
>>> wtay Tue May 8 20:45:42 2012
>>> [0]PETSC ERROR: Libraries linked from
>>> /home/wtay/Lib/petsc-3.2-dev_shared_debug/lib
>>> [0]PETSC ERROR: Configure run at Tue May 8 10:47:59 2012
>>> [0]PETSC ERROR: Configure options
>>> --with-mpi-dir=/opt/openmpi-1.5.3/
>>> --with-blas-lapack-dir=/opt/intelcpro-11.1.059/mkl/lib/em64t/
>>> --with-debugging=1 --download-hypre=1
>>> --prefix=/home/wtay/Lib/petsc-3.2-dev_shared_debug
>>> --known-mpi-shared=1 --with-shared-libraries
>>> [0]PETSC ERROR:
>>> ------------------------------------------------------------------------
>>> [0]PETSC ERROR: User provided function() line 0 in
>>> unknown directory unknown file
>>> --------------------------------------------------------------------------
>>> MPI_ABORT was invoked on rank 0 in communicator
>>> MPI_COMM_WORLD
>>> with errorcode 59.
>>>
>>> NOTE: invoking MPI_ABORT causes Open MPI to kill all MPI
>>> processes.
>>> You may or may not see output from other processes,
>>> depending on
>>> exactly when Open MPI kills them.
>>>
>>> Yours sincerely,
>>>
>>> TAY wee-beng
>>>
>>>
>>> On 5/5/2012 1:43 AM, Matthew Knepley wrote:
>>>> On Fri, May 4, 2012 at 5:42 PM, TAY wee-beng
>>>> <zonexo at gmail.com <mailto:zonexo at gmail.com>> wrote:
>>>>
>>>> Hi,
>>>>
>>>> I wonder if you people are interested to include my
>>>> ex29 fortran version in the petsc examples, which
>>>> can help people who are using fortran.
>>>>
>>>>
>>>> Yes, we will definitely include it. Please send the
>>>> source, and a representative output with run options.
>>>>
>>>> Thanks,
>>>>
>>>> Matt
>>>>
>>>> Thanks.
>>>>
>>>> Yours sincerely,
>>>>
>>>> TAY wee-beng
>>>>
>>>>
>>>> On 4/5/2012 9:28 PM, Matthew Knepley wrote:
>>>>> On Fri, May 4, 2012 at 3:24 PM, TAY wee-beng
>>>>> <zonexo at gmail.com <mailto:zonexo at gmail.com>> wrote:
>>>>>
>>>>>
>>>>> On 4/5/2012 9:16 PM, Barry Smith wrote:
>>>>>
>>>>> Do an hg pull and then run make in
>>>>> src/mat/interface/ftn-custom/
>>>>>
>>>>> Then it should link.
>>>>>
>>>>> Barry
>>>>>
>>>>> There was a E missing from the all caps
>>>>> name of the function.
>>>>>
>>>>> After hg pull, I did:
>>>>>
>>>>> cd src//mat/interface/ftn-custom/
>>>>>
>>>>> User at User-PC
>>>>> /cygdrive/c/temp/petsc-dev/src/mat/interface/ftn-custom
>>>>> $ make
>>>>> make[1]: Warning: File
>>>>> `/cygdrive/c/temp/petsc-dev/petsc-3.2-dev_win32_vs2008/lib/libpetsc.lib(zmatregf.o)'
>>>>> has modification time 787 s in the future
>>>>> make[1]: Nothing to be done for `libc'.
>>>>> make[1]: warning: Clock skew detected. Your
>>>>> build may be incomplete.
>>>>>
>>>>> But it still can't work.
>>>>>
>>>>>
>>>>> Something is messed up with the clock on this machine.
>>>>>
>>>>> HOWEVER, development requires certain basic skills
>>>>> in order to debug your work. We
>>>>> cannot be the ones debugging your code. Now
>>>>>
>>>>> nm $PETSC_ARCH/lib/libpetsc.a | grep -i
>>>>> MatNullSpaceRemove
>>>>>
>>>>> will look for the symbol.
>>>>>
>>>>> Matt
>>>>>
>>>>>
>>>>>
>>>>> On May 4, 2012, at 2:11 PM, Matthew
>>>>> Knepley wrote:
>>>>>
>>>>> On Fri, May 4, 2012 at 3:01 PM, TAY
>>>>> wee-beng<zonexo at gmail.com
>>>>> <mailto:zonexo at gmail.com>> wrote:
>>>>>
>>>>> On 4/5/2012 5:17 PM, Matthew Knepley
>>>>> wrote:
>>>>>
>>>>> On Fri, May 4, 2012 at 11:05 AM,
>>>>> TAY wee-beng<zonexo at gmail.com
>>>>> <mailto:zonexo at gmail.com>> wrote:
>>>>>
>>>>> On 4/5/2012 3:05 PM, Matthew
>>>>> Knepley wrote:
>>>>>
>>>>> On Fri, May 4, 2012 at 8:59
>>>>> AM, TAY
>>>>> wee-beng<zonexo at gmail.com
>>>>> <mailto:zonexo at gmail.com>> wrote:
>>>>>
>>>>> Hi,
>>>>>
>>>>> Is there anything else I can
>>>>> try to get it working right?
>>>>>
>>>>> The MatGetNullSpaceRemove() is
>>>>> missing.
>>>>>
>>>>> Where should I add
>>>>> MatGetNullSpaceRemove and what are
>>>>> its syntax? I googled but there's
>>>>> no results.
>>>>>
>>>>> Fixed in p;etsc-dev:
>>>>>
>>>>> http://www.mcs.anl.gov/petsc/petsc-dev/docs/manualpages/Mat/MatNullSpaceRemove.html
>>>>>
>>>>> I just compiled the updated petsc-dev
>>>>> but I got the same error msg when I use:
>>>>>
>>>>> call
>>>>> MatNullSpaceRemove(nullspace,b,PETSC_NULL,ierr)
>>>>>
>>>>> error LNK2019: unresolved external
>>>>> symbol MATNULLSPACEREMOVE referenced
>>>>> in function COMPUTERHS
>>>>> 1>c:\obj_tmp\ex29f\Debug\ex29f.exe :
>>>>> fatal error LNK1120: 1 unresolved
>>>>> externals
>>>>>
>>>>> That function is in:
>>>>>
>>>>> src/mat/interface/ftn-custom/zmatrixf.c
>>>>>
>>>>> Did that compile? Can you see the
>>>>> symbol in libpetsc.a?
>>>>>
>>>>> Matt
>>>>>
>>>>> Matt
>>>>>
>>>>> Thanks.
>>>>>
>>>>> Matt
>>>>>
>>>>> Thanks！
>>>>>
>>>>>
>>>>> On 2/5/2012 10:11 PM, Matthew
>>>>> Knepley wrote:
>>>>>
>>>>> On Wed, May 2, 2012 at
>>>>> 1:55 PM, TAY
>>>>> wee-beng<zonexo at gmail.com
>>>>> <mailto:zonexo at gmail.com>>
>>>>> wrote:
>>>>> Hi,
>>>>>
>>>>> I did a MatView and
>>>>> VecView on both C and
>>>>> Fortran, right after Mat
>>>>> and Vec assembly. I have
>>>>> attached the printout
>>>>> below. They are exactly
>>>>> the same, but yet the
>>>>> result is different in
>>>>> Neumann condition.
>>>>> However, the dirichlet
>>>>> condition gives the
>>>>> correct ans. Is there
>>>>> anything else that could
>>>>> be wrong even if the Mat
>>>>> and Vec are the same?
>>>>>
>>>>> Did you set the null space
>>>>> for the matrix when you
>>>>> have Neumann conditions?
>>>>>
>>>>> Yes, for the matrix, I set as:
>>>>>
>>>>> call
>>>>> MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,PETSC_NULL_OBJECT,nullspace,ierr)
>>>>>
>>>>> call
>>>>> MatSetNullSpace(jac,nullspace,ierr)
>>>>>
>>>>> call
>>>>> MatNullSpaceDestroy(nullspace,ierr)
>>>>>
>>>>> for the Vec,
>>>>>
>>>>> call
>>>>> MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,PETSC_NULL_OBJECT,nullspace,ierr)
>>>>>
>>>>> !call
>>>>> MatNullSpaceRemove(nullspace,b,PETSC_NULL,ierr)
>>>>>
>>>>> call
>>>>> MatNullSpaceDestroy(nullspace,ierr)
>>>>>
>>>>> MatNullSpaceRemove was comment
>>>>> out because there's error
>>>>> during linking
>>>>>
>>>>>
>>>>> Matt
>>>>>
>>>>> Thanks!
>>>>>
>>>>> Fortran:
>>>>>
>>>>> Matrix Object: 1 MPI processes
>>>>> type: seqaij
>>>>> row 0: (0, 2) (1, -1)
>>>>> (3, -1)
>>>>> row 1: (0, -1) (1, 3)
>>>>> (2, -1) (4, -1)
>>>>> row 2: (1, -1) (2, 2)
>>>>> (5, -1)
>>>>> row 3: (0, -1) (3, 3)
>>>>> (4, -1) (6, -1)
>>>>> row 4: (1, -1) (3, -1)
>>>>> (4, 4) (5, -1) (7, -1)
>>>>> row 5: (2, -1) (4, -1)
>>>>> (5, 3) (8, -1)
>>>>> row 6: (3, -1) (6, 2)
>>>>> (7, -1)
>>>>> row 7: (4, -1) (6, -1)
>>>>> (7, 3) (8, -1)
>>>>> row 8: (5, -1) (7, -1)
>>>>> (8, 2)
>>>>> Vector
>>>>> Object:Vec_0000000084000000_0
>>>>> 1 MPI processes
>>>>> type: mpi
>>>>> Process [0]
>>>>> 0.25
>>>>> 0.0205213
>>>>> 1.135e-005
>>>>> 0.0205213
>>>>> 0.00168449
>>>>> 9.31663e-007
>>>>> 1.135e-005
>>>>> 9.31663e-007
>>>>> 5.15289e-010
>>>>> Vector
>>>>> Object:Vec_0000000084000000_1
>>>>> 1 MPI processes
>>>>> type: mpi
>>>>> Process [0]
>>>>> 0.14924
>>>>> 0.0242397
>>>>> -0.0260347
>>>>> 0.0242397
>>>>> -0.0256192
>>>>> -0.0400102
>>>>> -0.0260347
>>>>> -0.0400102
>>>>> -0.0400102
>>>>> Press any key to continue
>>>>> . . .
>>>>>
>>>>> C:
>>>>>
>>>>> Matrix Object: 1 MPI processes
>>>>> type: seqaij
>>>>> row 0: (0, 2) (1, -1)
>>>>> (3, -1)
>>>>> row 1: (0, -1) (1, 3)
>>>>> (2, -1) (4, -1)
>>>>> row 2: (1, -1) (2, 2)
>>>>> (5, -1)
>>>>> row 3: (0, -1) (3, 3)
>>>>> (4, -1) (6, -1)
>>>>> row 4: (1, -1) (3, -1)
>>>>> (4, 4) (5, -1) (7, -1)
>>>>> row 5: (2, -1) (4, -1)
>>>>> (5, 3) (8, -1)
>>>>> row 6: (3, -1) (6, 2)
>>>>> (7, -1)
>>>>> row 7: (4, -1) (6, -1)
>>>>> (7, 3) (8, -1)
>>>>> row 8: (5, -1) (7, -1)
>>>>> (8, 2)
>>>>> Vector
>>>>> Object:Vec_0x1d3b000_0 1
>>>>> MPI processes
>>>>> type: mpi
>>>>> Process [0]
>>>>> 0.25
>>>>> 0.0205212
>>>>> 1.135e-05
>>>>> 0.0205212
>>>>> 0.00168449
>>>>> 9.31663e-07
>>>>> 1.135e-05
>>>>> 9.31663e-07
>>>>> 5.15288e-10
>>>>> Vector
>>>>> Object:Vec_0x1d3b000_1 1
>>>>> MPI processes
>>>>> type: mpi
>>>>> Process [0]
>>>>> 0.139311
>>>>> 0.0305751
>>>>> -0.0220633
>>>>> 0.0305751
>>>>> -0.0135158
>>>>> -0.042185
>>>>> -0.0220633
>>>>> -0.042185
>>>>> -0.058449
>>>>>
>>>>>
>>>>>
>>>>> Yours sincerely,
>>>>>
>>>>> TAY wee-beng
>>>>>
>>>>>
>>>>> On 1/5/2012 11:54 PM,
>>>>> Matthew Knepley wrote:
>>>>>
>>>>> On Tue, May 1, 2012 at
>>>>> 5:48 PM, TAY
>>>>> wee-beng<zonexo at gmail.com
>>>>> <mailto:zonexo at gmail.com>>
>>>>> wrote:
>>>>> 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.
>>>>>
>>>>> This is very simple.
>>>>> You have an error in
>>>>> your code. Checking it
>>>>> is very simple: run
>>>>> the code and
>>>>> break in
>>>>> MatSetValues(). Make
>>>>> sure ex29 makes calls
>>>>> with exactly the same
>>>>> indices as your ex29f.
>>>>>
>>>>> Matt
>>>>>
>>>>> 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
>>>>>
>>>>>
>>>>>
>>>>> --
>>>>> 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/20120603/5ccf914c/attachment-0001.html>
-------------- next part --------------
! Concepts: KSP^solving a system of linear equations
! Concepts: KSP^Laplacian, 2d
! Processors: n
!Added by Wee Beng Tay (W.B.Tay at tudelft.nl).
!conversion from ex29 C version to Fortran
!Inhomogeneous Laplacian in 2D. Modeled by the partial differential equation
! -div \rho grad u = f, 0 < x,y < 1,
!with forcing function
! f = e^{-x^2/\nu} e^{-y^2/\nu}
!with Dirichlet boundary conditions
! u = f(x,y) for x = 0, x = 1, y = 0, y = 1
!or pure Neumman boundary conditions
!This uses multigrid to solve the linear system
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/petscsys.h>
!#include <finclude/petscvec.h>
!#include <finclude/petscmat.h>
!#include <finclude/petscpc.h>
!#include <finclude/petscksp.h>
!#include <finclude/petscdmda.h>
!#include <finclude/petscdmda.h90>
#include "finclude/petsc.h90"
PetscErrorCode ierr
DM da
KSP ksp
Vec x
external ComputeRHS,ComputeMatrix
PetscInt i1,i3,bc_cond !>bc_cond = boundary condition = 1/2 = dirichlet/neumann
!>have to change bc_cond at 3 locations
call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
bc_cond = 1
i3 = -3
i1 = 1
call KSPCreate(PETSC_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/petscsys.h>
!#include <finclude/petscvec.h>
!#include <finclude/petscmat.h>
!#include <finclude/petscpc.h>
!#include <finclude/petscksp.h>
!#include <finclude/petscdmda.h90>
!#include <finclude/petscdmda.h>
#include "finclude/petsc.h90"
PetscErrorCode ierr
PetscInt i,j,ij,mx,my,xm,ym,xs,ys,bc_cond,status,myid,ij_sta,ij_end
PetscScalar h,nu,rho
PetscScalar Hx,Hy
PetscScalar,pointer :: array(:,:) !Problem with using DMDAVecGetArrayF90 and DMDAVecRestoreArrayF90 below in
Vec x,b !linux with ifort, so VecSetValues was used. Worked in VS2008 though
DM da
MatNullSpace nullspace !>neumann BC
integer, allocatable:: ij_array(:)
real(8), allocatable :: arrays(:)
nu = 0.1
bc_cond = 1
call MPI_COMM_RANK(MPI_COMM_WORLD, myid, ierr)
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 VecGetOwnershipRange(b,ij_sta,ij_end,ierr)
allocate (ij_array(ij_sta:ij_end-1), STAT=status)
if (status/=0) STOP "Cannot allocate memory 1"
allocate (arrays(ij_sta:ij_end-1), STAT=status)
if (status/=0) STOP "Cannot allocate memory 2"
!Section below not usable in linux ifort yet - only works in VS2008 + ifort
!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)
!Section above not usable in linux ifort yet - only works in VS2008 + ifort
ij = ij_sta - 1
do j = ys,ys+ym-1
do i = xs,xs+xm-1
ij = ij + 1
ij_array(ij) = ij
arrays(ij) = exp(-(i*Hx)*(i*Hx)/nu)*exp(-(j*Hy)*(j*Hy)/nu)*Hx*Hy
end do
end do
call VecSetValues(b,ij_end-ij_sta,ij_array,arrays,INSERT_VALUES,ierr)
call VecAssemblyBegin(b,ierr)
call VecAssemblyEnd(b,ierr)
!call VecView(b,PETSC_VIEWER_STDOUT_WORLD,ierr)
if (bc_cond == 2) then
call MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,PETSC_NULL_OBJECT,nullspace,ierr)
call MatNullSpaceRemove(nullspace,b,PETSC_NULL_OBJECT,ierr)
call MatNullSpaceDestroy(nullspace,ierr)
end if
end subroutine ComputeRHS
subroutine ComputeMatrix(da,x,JJ,jac,str,ierr)
implicit none
!#include <finclude/petscsys.h>
!#include <finclude/petscvec.h>
!#include <finclude/petscmat.h>
!#include <finclude/petscdmda.h>
!#include <finclude/petscdmda.h90>
#include "finclude/petsc.h90"
Mat jac,JJ
PetscErrorCode ierr
DM da
PetscInt i,j,k,mx,my,xm,num
PetscInt ym,xs,ys,i1,i5,bc_cond
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
bc_cond = 1
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==0 .or. j==0 .or. i==mx-1 .or. j==my-1) then
if (bc_cond == 1) then
v(1) = 2.0*rho*(HxdHy + HydHx)
call MatSetValuesStencil(jac,i1,row,i1,row,v,INSERT_VALUES,ierr)
else if (bc_cond == 2) then
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)
col(MatStencil_i,num) = i
col(MatStencil_j,num) = j
!num = num + 1
call MatSetValuesStencil(jac,i1,row,num,col,v,INSERT_VALUES,ierr)
end if
else
v(1) = -rho*HxdHy
col(MatStencil_i,1) = i
col(MatStencil_j,1) = j-1
v(2) = -rho*HydHx
col(MatStencil_i,2) = i-1
col(MatStencil_j,2) = j
v(3) = 2.0*rho*(HxdHy + HydHx)
col(MatStencil_i,3) = i
col(MatStencil_j,3) = j
v(4) = -rho*HydHx
col(MatStencil_i,4) = i+1
col(MatStencil_j,4) = j
v(5) = -rho*HxdHy
col(MatStencil_i,5) = i
col(MatStencil_j,5) = j+1
call MatSetValuesStencil(jac,i1,row,i5,col,v,INSERT_VALUES,ierr)
end if
end do
end do
call MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr)
call MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr)
!call MatView(jac,PETSC_VIEWER_STDOUT_WORLD,ierr)
if (bc_cond == 2) then
call MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,PETSC_NULL_OBJECT,nullspace,ierr)
call MatSetNullSpace(jac,nullspace,ierr)
call MatNullSpaceDestroy(nullspace,ierr)
end if
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
More information about the petsc-users
mailing list