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
