[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