Re: [petsc-users] accessing DMDA Vec ghost values
Thanks Matt and Dave for the explanation, that is very helpful.
Best
Sean
On Thu, May 12, 2016 at 8:58 AM -0700, "Dave May"
<dave.mayhe...@gmail.com<mailto:dave.mayhe...@gmail.com>> wrote:
Matt beat me to the punch... :D
Anyway, here is my more detailed answer.
Thanks! Somehow I missed DM{Get,Create}LocalVector(). BTW what is the
difference between the Get and Create versions? It is not obvious from the
documentation.
The DMDA contains a pool of vectors (both local and global) which can be
re-used by the user. This avoids the need to continually allocate and
deallocate memory. Thus, the Get methods are simply an optimization.
The Get methods retrieve from the pool, a vector which isn't currently in use.
In this case, You can think of Restore as returning the vector back to the pool
to be used somewhere else.
If all vectors in the pool are in use, a new one will be allocated for you. In
this case, Restore will actually deallocate memory.
Since Get methods may return vectors which have been used else where in the
code, you should always call VecZeroEntries() on them.
The Create methods ALWAYS allocate new memory and thus you ALWAYS need to call
Destroy on them. Vectors obtained via VecCreate() will always be initialized
with 0's.
Thanks,
Dave
Also, can you explain the difference between DMDAVecGetArrayDOF and
DMDAVecGetArrayDOFRead?
Thanks again,
Sean
Thanks,
Dave
Thanks very much!
Sean Dettrick
Re: [petsc-users] accessing DMDA Vec ghost values
From: <petsc-users-boun...@mcs.anl.gov<mailto:petsc-users-boun...@mcs.anl.gov>>
on behalf of Dave May <dave.mayhe...@gmail.com<mailto:dave.mayhe...@gmail.com>>
Date: Thursday, May 12, 2016 at 2:48 AM
To: Sean Dettrick
<sdettr...@trialphaenergy.com<mailto:sdettr...@trialphaenergy.com>>
Cc: "petsc-users@mcs.anl.gov<mailto:petsc-users@mcs.anl.gov>"
<petsc-users@mcs.anl.gov<mailto:petsc-users@mcs.anl.gov>>
Subject: Re: [petsc-users] accessing DMDA Vec ghost values
On 12 May 2016 at 10:42, Sean Dettrick
<sdettr...@trialphaenergy.com<mailto:sdettr...@trialphaenergy.com>> wrote:
Hi,
When discussing DMDAVecGetArrayDOF etc in section 2.4.4, the PETSc 3.7 manual
says "The array is accessed using the usual global indexing on the entire grid,
but the user may only refer to the local and ghost entries of this array as all
other entries are undefined”.
OK so far. But how to access the ghost entries?
With a 2D DMDA, I can do this OK:
PetscIntxs,xm,ys,ym;
ierr=DMDAGetCorners(da,,,0,,,0);CHKERRQ(ierr);
PetscScalar ***es;
ierr=DMDAVecGetArrayDOF(da,Es,);CHKERRQ(ierr);
for (int j=ys; j < ys+ym; j++) {
for (int i=xs; i < xs+xm;i++) {
es[j][i][0]=1.;
es[j][i][1]=1.;
}
}
ierr=DMDAVecRestoreArrayDOF(da,Es,);CHKERRQ(ierr);
But if I replace DMDAGetCorners with DMDAGetGhostCorners, then the code crashes
with a seg fault, presumably due to out of bounds memory access.
Is that supposed to happen?
If you created the vector Es using the function DM{Get,Create}GlobalVector(),
then the answer is yes.
What’s the remedy?
If you want to access the ghost entries, you need to create the vector using
the function DM{Get,Create}LocalVector().
Thanks! Somehow I missed DM{Get,Create}LocalVector(). BTW what is the
difference between the Get and Create versions? It is not obvious from the
documentation.
Also, can you explain the difference between DMDAVecGetArrayDOF and
DMDAVecGetArrayDOFRead?
Thanks again,
Sean
Thanks,
Dave
Thanks very much!
Sean Dettrick
[petsc-users] accessing DMDA Vec ghost values
Hi,
When discussing DMDAVecGetArrayDOF etc in section 2.4.4, the PETSc 3.7 manual
says "The array is accessed using the usual global indexing on the entire grid,
but the user may only refer to the local and ghost entries of this array as all
other entries are undefined”.
OK so far. But how to access the ghost entries?
With a 2D DMDA, I can do this OK:
PetscIntxs,xm,ys,ym;
ierr=DMDAGetCorners(da,,,0,,,0);CHKERRQ(ierr);
PetscScalar ***es;
ierr=DMDAVecGetArrayDOF(da,Es,);CHKERRQ(ierr);
for (int j=ys; j < ys+ym; j++) {
for (int i=xs; i < xs+xm;i++) {
es[j][i][0]=1.;
es[j][i][1]=1.;
}
}
ierr=DMDAVecRestoreArrayDOF(da,Es,);CHKERRQ(ierr);
But if I replace DMDAGetCorners with DMDAGetGhostCorners, then the code crashes
with a seg fault, presumably due to out of bounds memory access.
Is that supposed to happen? What’s the remedy?
Thanks very much!
Sean Dettrick
Re: [petsc-users] Petsc Draw to pixmap
No, it is my fault, when I was unable to download from afterstep.org<http://afterstep.org> earlier I didn’t realize that it was probably my work firewall blocking me. Thus sending me on a wild goose chase with macports and home-brew and github version, which didn’t compile. But now I’ve tried to follow your advice, and have a new problem. I installed Afterimage from your suggested site, with suggested X lib flags: ./configure --x-includes=/opt/X11/include -x-libraries=/opt/X11/lib make make install No errors. Then re-configured and built petsc: export PETSC_DIR=/Users/sdettrick/libs/petsc-3.6.1 export PETSC_ARCH=mac-gcc-hdf5-mpich-afterimage ./configure --with-afterimage --download-hdf5 --download-mpich --PETSC_ARCH=mac-gcc-hdf5-mpich-afterimage --PETSC_DIR=/Users/sdettrick/libs/petsc-3.6.1 make PETSC_DIR=/Users/sdettrick/libs/petsc-3.6.1 PETSC_ARCH=mac-gcc-hdf5-mpich-afterimage all No errors. Then try to run an example with image output: cd src/dm/examples/tutorials/ make ex5 ./ex5 ./ex5 -draw_save test.png The example generates the error you mentioned, which it didn’t before installing afterimage: X Error of failed request: BadValue (integer parameter out of range for operation) Major opcode of failed request: 1 (X_CreateWindow) Value in failed request: 0x40 Serial number of failed request: 7 Current serial number in output stream: 14 Browsing through configure.log, it seems that petsc is using the same X11 library. Did I miss something? Thanks! Sean On Dec 2, 2015, at 11:10 PM, Barry Smith <bsm...@mcs.anl.gov<mailto:bsm...@mcs.anl.gov>> wrote: I am such an idiot. No I did build it from source I forgot the comment in afterimage.py After image is available from http://www.afterstep.org/afterimage/getcode.php # # It is used by the PetscDrawSetSave() routine to save X windows graphics to files # # If installing on an Apple make sure to read the details on PetscDrawSetSave manual # page before installing # from that page I found If X windows generates an error message about X_CreateWindow() failing then Afterimage was installed without X windows. Reinstall Afterimage using the ./configure flags --x-includes=/pathtoXincludes --x-libraries=/pathtoXlibraries For example under Mac OS X Mountain Lion --x-includes=/opt/X11/include -x-libraries=/opt/X11/lib Did you try all this stuff? Barry On Dec 3, 2015, at 1:00 AM, Sean Dettrick <sdettr...@trialphaenergy.com<mailto:sdettr...@trialphaenergy.com>> wrote: Hi Barry, Thanks for the amazingly fast reply. I installed home-brew but then was unable to find afterstep available as a package. The closest thing I could find was asterm (after step terminal emulator. Do you remember if you installed it from a canonical repo? I did find afterstep on github, but had trouble installing it on the mac. Meanwhile, a colleague has given me a python script, so I’ll probably use that instead for now. Thanks again, Sean On Dec 2, 2015, at 1:45 PM, Barry Smith <bsm...@mcs.anl.gov<mailto:bsm...@mcs.anl.gov>> wrote: Sean, I have never been able to build from that source either. I used homebrew to install it and its worked for several years. I also tried to see if one could pull the image from the X pixel map and it looked like a huge project (the pixmap stuff is not as trivial as one would think it would be) so I gave up on that and started using afterimage (which essentially does all the for you). Barry On Dec 2, 2015, at 3:39 PM, Sean Dettrick <sdettr...@trialphaenergy.com<mailto:sdettr...@trialphaenergy.com>> wrote: Hi, I’d like to have petsc draw output to a pixmap, but have been unable to install afterstep on my Mac. Direct download of the tar files from http://www.afterstep.org/afterimage/getcode.php failed, as did cvs checkout. It seems the target sites may not exist at present. Is there another way to install afterstep? Alternatively, is there a way to manually redirect the X11 output using something like XCreatePixmap? Thanks! Sean Dettrick
Re: [petsc-users] Petsc Draw to pixmap
On Dec 3, 2015, at 5:01 AM, Barry Smith <bsm...@mcs.anl.gov<mailto:bsm...@mcs.anl.gov>> wrote: On Dec 3, 2015, at 5:54 AM, Matthew Knepley <knep...@gmail.com<mailto:knep...@gmail.com>> wrote: On Thu, Dec 3, 2015 at 2:06 AM, Sean Dettrick <sdettr...@trialphaenergy.com<mailto:sdettr...@trialphaenergy.com>> wrote: No, it is my fault, when I was unable to download from afterstep.org<http://afterstep.org> earlier I didn’t realize that it was probably my work firewall blocking me. Thus sending me on a wild goose chase with macports and home-brew and github version, which didn’t compile. But now I’ve tried to follow your advice, and have a new problem. I installed Afterimage from your suggested site, with suggested X lib flags: ./configure --x-includes=/opt/X11/include -x-libraries=/opt/X11/lib Did you check that it worked? Configure can (semi-)silently turn things off. Yes, I think Matt is right. This is the error message you get if afterimage was installed without X. Check through all the logs and build stuff of afterimage. Thanks! The thing finally works. I had to clean everything and do fresh installs, but now it is fine. Summarizing my experience of afterimage on a mac in case anybody else wants to do it: Afterimage on the mac is in theory available via macports as part of the afterstep package, but it failed to build on my Mac (Yosemite). It is not available in home-brew. The current github version of afterimage failed to build on my mac. The version from http://www.afterstep.org/afterimage/getcode.php DOES work, configured as: ./configure --x-includes=/opt/X11/include --x-libraries=/opt/X11/lib --with-x make make install Then petsc is configured as: ./configure --with-afterimage --download-hdf5 --download-mpich The example src/dm/examples/tutorials/ex5.c then works: ./ex5 -draw_save Cheers, Sean Barry Thanks, Matt make make install No errors. Then re-configured and built petsc: export PETSC_DIR=/Users/sdettrick/libs/petsc-3.6.1 export PETSC_ARCH=mac-gcc-hdf5-mpich-afterimage ./configure --with-afterimage --download-hdf5 --download-mpich --PETSC_ARCH=mac-gcc-hdf5-mpich-afterimage --PETSC_DIR=/Users/sdettrick/libs/petsc-3.6.1 make PETSC_DIR=/Users/sdettrick/libs/petsc-3.6.1 PETSC_ARCH=mac-gcc-hdf5-mpich-afterimage all No errors. Then try to run an example with image output: cd src/dm/examples/tutorials/ make ex5 ./ex5 ./ex5 -draw_save test.png The example generates the error you mentioned, which it didn’t before installing afterimage: X Error of failed request: BadValue (integer parameter out of range for operation) Major opcode of failed request: 1 (X_CreateWindow) Value in failed request: 0x40 Serial number of failed request: 7 Current serial number in output stream: 14 Browsing through configure.log, it seems that petsc is using the same X11 library. Did I miss something? Thanks! Sean On Dec 2, 2015, at 11:10 PM, Barry Smith <bsm...@mcs.anl.gov<mailto:bsm...@mcs.anl.gov>> wrote: I am such an idiot. No I did build it from source I forgot the comment in afterimage.py After image is available from http://www.afterstep.org/afterimage/getcode.php # # It is used by the PetscDrawSetSave() routine to save X windows graphics to files # # If installing on an Apple make sure to read the details on PetscDrawSetSave manual # page before installing # from that page I found If X windows generates an error message about X_CreateWindow() failing then Afterimage was installed without X windows. Reinstall Afterimage using the ./configure flags --x-includes=/pathtoXincludes --x-libraries=/pathtoXlibraries For example under Mac OS X Mountain Lion --x-includes=/opt/X11/include -x-libraries=/opt/X11/lib Did you try all this stuff? Barry On Dec 3, 2015, at 1:00 AM, Sean Dettrick <sdettr...@trialphaenergy.com<mailto:sdettr...@trialphaenergy.com>> wrote: Hi Barry, Thanks for the amazingly fast reply. I installed home-brew but then was unable to find afterstep available as a package. The closest thing I could find was asterm (after step terminal emulator. Do you remember if you installed it from a canonical repo? I did find afterstep on github, but had trouble installing it on the mac. Meanwhile, a colleague has given me a python script, so I’ll probably use that instead for now. Thanks again, Sean On Dec 2, 2015, at 1:45 PM, Barry Smith <bsm...@mcs.anl.gov<mailto:bsm...@mcs.anl.gov>> wrote: Sean, I have never been able to build from that source either. I used homebrew to install it and its worked for several years. I also tried to see if one could pull the image from the X pixel map and it looked like a huge project (the pixmap stuff is not as trivial as one would think it would be) so I gave up on that and started using afterimage (which essentially does all the for you). Barry On Dec 2
[petsc-users] Petsc Draw to pixmap
Hi, I’d like to have petsc draw output to a pixmap, but have been unable to install afterstep on my Mac. Direct download of the tar files from http://www.afterstep.org/afterimage/getcode.php failed, as did cvs checkout. It seems the target sites may not exist at present. Is there another way to install afterstep? Alternatively, is there a way to manually redirect the X11 output using something like XCreatePixmap? Thanks! Sean Dettrick
Re: [petsc-users] Petsc Draw to pixmap
Hi Barry, Thanks for the amazingly fast reply. I installed home-brew but then was unable to find afterstep available as a package. The closest thing I could find was asterm (after step terminal emulator. Do you remember if you installed it from a canonical repo? I did find afterstep on github, but had trouble installing it on the mac. Meanwhile, a colleague has given me a python script, so I’ll probably use that instead for now. Thanks again, Sean > On Dec 2, 2015, at 1:45 PM, Barry Smith <bsm...@mcs.anl.gov> wrote: > > > Sean, > > I have never been able to build from that source either. > > I used homebrew to install it and its worked for several years. > > I also tried to see if one could pull the image from the X pixel map and it > looked like a huge project (the pixmap stuff is not as trivial as one would > think it would be) so I gave up on that and started using afterimage (which > essentially does all the for you). > > Barry > >> On Dec 2, 2015, at 3:39 PM, Sean Dettrick <sdettr...@trialphaenergy.com> >> wrote: >> >> Hi, >> >> I’d like to have petsc draw output to a pixmap, but have been unable to >> install afterstep on my Mac. Direct download of the tar files from >> http://www.afterstep.org/afterimage/getcode.php failed, as did cvs checkout. >> It seems the target sites may not exist at present. >> >> Is there another way to install afterstep? >> >> Alternatively, is there a way to manually redirect the X11 output using >> something like XCreatePixmap? >> >> Thanks! >> Sean Dettrick >
[petsc-users] PETSc on Blue Gene/Q ?
Dear Petsc Users, Does anyone have experience with running a PETSc app or benchmark on Blue Gene/Q? Does PETSc compile and run? Does it link to the IBM MASS libraries? How does a single Blue Gene CPU compare to a Xeon? Does your PETSc code scale well on BGQ? Thanks! Sean Dettrick
[petsc-users] Problems with the intel compilers
Intel optimization is slow, and is turned on by default. If you turn off intel optimization completely with -O0, then you'll probably find compilation speeds up to be similar to gnu. Is the increased compilation time of intel worth it? We find intel programs are usually faster than gnu ones. And the intel debugger idb is in many ways more powerful than gdb. But on the other hand we have a fortran based CG solver where gfortran compiled code is much faster than the intel one. Cheers, Sean On Oct 27, 2009, at 3:54 PM, Dominik Szczerba wrote: Thanks Ando, I know of the slight speed advantage in some cases on the Intel side, but the question was more why the compilation takes so much longer? Dominik Andreas Grassl wrote: Hi Dominik, Dominik Szczerba wrote: I have noticed already a longer while ago: setting up PETSc with the Intel compilers (--with-cc=icc --with-cxx=icpc) takes an order of magnitude longer than with the native GNU. The step taking most of the time is configuring mpich. I have tried to configure mpich separately and indeed, with gnu it is a breeze, with intel 10x slower. In both cases, linux both 32 and 64 bit, Ubuntu 9.04 and debian testing. Intel compilers 10.x and 11.x. I just wanted to ask opinion if anybody has similar observations and/or finds using intel worthwhile at all. On itanium machines intel compilers produce much faster code compared to gnu-compilers. On other machines I did not notice a big speed difference, but I think intel compilers produce slightly faster code as well. Cheers, ando
Performance degradation after upgrade to 3.0.0
This problem occurs when the MS user sends an attachment to an email which has rich text formatting. Format the email in plain text, attach, and try again. S On Feb 2, 2009, at 5:30 AM, Barry Smith wrote: There may be some option (well hidden) in outlook web access that lets you deselect winmail.dat but it is using it by default. Barry On Feb 2, 2009, at 3:59 AM, Billy Ara?jo wrote: I didn't send winmail.dat. Maybe because was sent from outlook web access, I don't know... :) Billy. -Original Message- From: petsc-users-bounces at mcs.anl.gov on behalf of Barry Smith Sent: Sun 2/1/2009 11:24 PM To: PETSc users list Subject: Re: Performance degradation after upgrade to 3.0.0 winmail.dat? Come on, be serious, we'd like to help you but not all of us are chained to Bill Gates nightstand. Barry On Feb 1, 2009, at 3:01 PM, Billy Ara?jo wrote: Hi, Here are the files that contain all the information below. You can see that the matrices are the same, unless I am missing something here. :) (I am resending them in compressed format). Regards, Billy. -Original Message- From: petsc-users-bounces at mcs.anl.gov on behalf of Matthew Knepley Sent: Sun 2/1/2009 6:07 PM To: PETSc users list Subject: Re: Performance degradation after upgrade to 3.0.0 In order to determine what is happening you first should: 1) Confirm that the solver setup is identical using -ksp_view for both versions 2) Determine that the matrices are identical. Output both matrices using MatView() with a PetscBinaryViewer. You can diff the files, and you can also solver both matrices using KSP ex10. 3) Look at the residuals using -ksp_monitor. If they are different, something else has changed. Matt On Sun, Feb 1, 2009 at 6:32 AM, Billy Ara?jo billy at dem.uminho.pt wrote: Hi, I have also verified that there has been a degradation of performance using the new 3.0 version: This is my function calling PETSc: KSP ksp; PC pc; KSPCreate (PETSC_COMM_WORLD, ksp); KSPSetOperators (ksp, *A, *A, DIFFERENT_NONZERO_PATTERN); KSPSetType (ksp, KSPFGMRES); KSPGetPC (ksp, pc); PCSetType (pc, PrecondProc); KSPSetInitialGuessNonzero (ksp, PETSC_TRUE); KSPSetTolerances (ksp, 1E-50, maxtol, PETSC_DEFAULT, maxiter); KSPSetFromOptions (ksp); KSPSolve (ksp, *b, *x); KSPGetIterationNumber (ksp, iter); KSPGetResidualNorm (ksp, res); KSPDestroy (ksp); with previous version 2.3.3-p6: Number of iterations: 42 Residual: +7.073781E-13 Time: +8.615024E-03 now: Number of iterations: 500 Residual: +2.746161E-05 Time: +1.026870E-01 It is reaching maximum number of iterations. The only thing I changed was: MatSetOption (*A, MAT_SYMMETRIC); to MatSetOption (*A, MAT_SYMMETRIC, PETSC_TRUE); I think I didn't change anything else. Regards, Billy. -- 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 winmail.dat winmail.dat
process id in DA
On Nov 12, 2008, at 4:24 PM, Shao-Ching Huang wrote:
Hi,
Is there a PETSc call to obtain the process coordinates in a process
grid, partitioned by a DA? (Something like what MPI_Cart_coords does.)
Thanks,
Shao-Ching
A while ago (May 2006) I found that the PETSc DA has its processes
ordered in the TRANSPOSE way to that of MPI Cartesian communicators,
so you might have to consider mapping back and forth between the two
cartesian representations.
Barry said that he might change the DA so that the two agree with each
other, but I don't know if he ever did...
Does anybody know the status of that?
In case it helps, below is part of our exchange on the matter.
Sean
---
Ok, so it is laid-out like a matrix, not like x, y coordinates.
Will put the change to DA on the list of things to do.
Barry
On Thu, 18 May 2006, Sean Dettrick wrote:
Barry Smith wrote:
Bill, does the MPI standard dictate this decomposition or
could different implementations do it the opposite way?
Then we'd have to make the DA logic a bit more complicated.
I don't have a copy of the standard, but to quote page 255 of MPI,
the complete reference by Snir et al:
Row-major numbering is always used for the processes in a Cartesian
structure.
Their diagram in figure 6.1 matches my code output for coords
couplets (i,j):
0 1 2 3
(0,0) (0,1) (0,2) (0,3)
4 5 6 7
(1,0) (1,1) (1,2) (1,3)
8 9 10 11
(2,0) (2,1) (2,2) (2,3)
By the way I agree with you, I *should* be able to swap the x and y
myself. Just haven't had much luck yet in that regard.
Sean
On Thu, 18 May 2006, Sean Dettrick wrote:
Barry Smith wrote:
On Thu, 18 May 2006, Sean Dettrick wrote:
Hi Barry,
the order is determined by MPI_Cart_create.
Do you mean that MPI_Cart_create() orders across the 2nd (y-axis)
fastest and then the first (x-axis)? Hmmm, maybe we should change
the
DA? Changing it once and for all (not supporting both) is probably
not a big deal and shouldn't break much (I hope).
Hi Barry,
it depends, what do you call x and what do you call y?
MPI_Cart_coords returns a vector, coords - I tend to say x is
coords[0], y is coords[1] and z is coords[2]. For what it's worth,
there's a short code appended to this email, which produces:
rank = 0 has Cartesian coords = { 0, 0 }
rank = 1 has Cartesian coords = { 0, 1 }
rank = 2 has Cartesian coords = { 1, 0 }
rank = 3 has Cartesian coords = { 1, 1 }
rank = 0 has DA range x=[0,50) and y=[0,50)
rank = 1 has DA range x=[50,100) and y=[0,50)
rank = 2 has DA range x=[0,50) and y=[50,100)
rank = 3 has DA range x=[50,100) and y=[50,100)
I don't completely understand what goes wrong. Is it because
YOUR
application orders the processors related to geometry in the
following way?
^ y direction
|
2 5 8
1 4 7
0 3 6
- x direction
Or is this something inherent in MPI_Cart_create?
For my interpretation of x and y, MPI_Cart_create produces the
above layout. But if I said x=coords[1] and y=coords[0], then it
would match the one below.
PETSc does it so
^ y direction
|
6 7 8
3 4 5
0 1 2
- x direction
Code and makefile attached ... hopefully within the message size
limit.
Just make cartcommtest.
Sean
mixed matrix type?
Hi, I have a sparse N*N matrix generated from a DA and a 5 point stencil, with a total of approx 5*N non-zero entries. Now I would like to extend this matrix by adding a smaller M*M dense matrix to the bottom right hand side, i.e. so that there is a dense square in the bottom right hand corner of the otherwise sparse matrix. The total number of new non-zero entries, M*M, is comparable to the total number of old entries, 5*N. On top of this, there would be a small number of non- zero entries in the new upper-right and lower-left rectangular portions of the matrix, due to coupling of the two systems. The new total matrix size (including zeroes) would be (N+M)*(N+M). Can anybody recommend a Mat type to store the new matrix? One possibility I was thinking of was to establish the original sparse Mat with a DA in a sub-communicator (with half the CPUs), and get the ownership range with MatGetOwnershipRange. Then in the petsc_comm_world communicator, the complete matrix could be constructed (by element-wise copying the old one I suppose), and the ownership range could be maintained manually. Does this sound like a reasonable strategy? I would very much appreciate any suggestions or advice. Thanks, Sean Dettrick
mixed matrix type?
Thanks Matt, thanks Barry, much appreciated. Best, Sean
How to efficiently change just the diagonal vector in a matrix at every time step
On May 10, 2008, at 10:38 AM, Ben Tay wrote: Hi Sean, Maybe for me, I can just insert vector diagonal D into the matrix A and call Assembly and KSP at every time step. Should that be better since there is no need to copy A into B? Thanks! Yes, that sounds like it would be faster. I suppose you would use INSERT_VALUES rather than ADD_VALUES. In my case I copy the whole matrix because constructing the time- constant part of the diagonal is very complicated. But now that you mention it, I could store both the time-constant and time-varying diagonal components in two separate Vecs, and only have one Mat, and then do MatDiagonalSet() twice at each timestep - the first time with INSERT_VALUES, the second with ADD_VALUES. That sounds like it would be faster. Thanks to you too, Sean On Fri, May 9, 2008 at 8:50 PM, Sean Dettrick sdettrick at gmail.com wrote: One way to do it is to have two Mats, A and B, and a Vec, D, to store the diagonal. A is constructed only on the first step. On subsequent steps, A is copied into B, and then D is added to the diagonal: ierr = MatCopy( A, B, SAME_NON_ZERO_PATTERN ); ierr = MatDiagonalSet( B, D, ADD_VALUES ); The KSP uses B as the matrix, not A. I don't know if this approach is efficient or not. Can anybody comment? Thanks, Sean On May 9, 2008, at 2:33 PM, Ben Tay wrote: Hi, I have a matrix and I inserted all the relevant values during the 1st step. I'll then solve it. For the subsequent steps, I only need to change the diagonal vector of the matrix before solving. I wonder how I can do it efficiently. Of course, the RHS vector also change but I've not included them here. I set these at the 1st step: call KSPSetOperators (ksp_semi_x,A_semi_x,A_semi_x,SAME_NONZERO_PATTERN,ierr) call KSPGetPC(ksp_semi_x,pc_semi_x,ierr) ksptype=KSPRICHARDSON call KSPSetType(ksp_semi_x,ksptype,ierr) ptype = PCILU call PCSetType(pc_semi_x,ptype,ierr) call KSPSetFromOptions(ksp_semi_x,ierr) call KSPSetInitialGuessNonzero(ksp_semi_x,PETSC_TRUE,ierr) tol=1.e-5 call KSPSetTolerances (ksp_semi_x ,tol ,PETSC_DEFAULT_DOUBLE_PRECISION ,PETSC_DEFAULT_DOUBLE_PRECISION,PETSC_DEFAULT_INTEGER,ierr) and what I did at the subsequent steps is: do II=1,total call MatSetValues(A_semi_x,1,II,1,II,new_value,INSERT_VALUES,ierr) end do call MatAssemblyBegin(A_semi_x,MAT_FINAL_ASSEMBLY,ierr) call MatAssemblyEnd(A_semi_x,MAT_FINAL_ASSEMBLY,ierr) call KSPSolve(ksp_semi_x,b_rhs_semi_x,xx_semi_x,ierr) I realise that the answers are slightly different as compared to calling all the options such as KSPSetType, KSPSetFromOptions, KSPSetTolerances at every time step. Should that be so? Is this the best way? Also, I can let the matrix be equal at every time step by fixing the delta_time. However, it may give stability problems. I wonder how expensive is these type of value changing and assembly for a matrix? Thank you very much. Regards. -- next part -- An HTML attachment was scrubbed... URL: http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20080512/13c5d22e/attachment.htm
DA question
To elaborate on Matt's suggestion, a staggered grid/Yee mesh code could use a single DA with one degree-of-freedom per component of H and E. The extra overlap required for staggered guard cells at the domain boundaries could be dealt with by having a bigger-than-usual stencil width. For the 2nd order 3D case, this suggests the DACreate3d routine would have arguments dof=6, s=2, and stencil_type=DA_STENCIL_STAR. It is just a suggestion - I have not tried it. Sean On Wed, Apr 9, 2008 at 5:06 PM, Amit.Itagi at seagate.com wrote: Randy, I guess, since you are doing a frequency domain calculation, you eventually end up with a single matrix equation. I am planning to work in the time domain. Will that change things ? Thanks Rgds, Amit Randall Mackie rlmackie862 at gmai l.com To Sent by: petsc-users at mcs.anl.gov owner-petsc-users cc @mcs.anl.gov No Phone Info Subject Available Re: DA question 04/09/2008 04:09 PM Please respond to petsc-users at mcs.a nl.gov Hi Amit, Why do you need two staggered grids? I do EM finite difference frequency domain modeling on a staggered grid using just one DA. Works perfectly fine. There are some grid points that are not used, but you just set them to zero and put a 1 on the diagonal of the coefficient matrix. Randy Amit.Itagi at seagate.com wrote: Hi Berend, A detailed explanation of the finite difference scheme is given here : http://en.wikipedia.org/wiki/Finite-difference_time-domain_method Thanks Rgds, Amit Berend van Wachem berend at chalmers. se To Sent by: petsc-users at mcs.anl.gov owner-petsc-users cc @mcs.anl.gov No Phone Info Subject Available Re: DA question 04/09/2008 02:59 PM Please respond to petsc-users at mcs.a nl.gov Dear Amit, Could you explain how the two grids are attached? I am using multiple DA's for multiple structured grids glued together. I've done the gluing with setting up various IS objects. From the multiple DA's, one global variable vector is formed. Is that what you are looking for? Best regards, Berend. Amit.Itagi at seagate.com wrote: Hi, Is it possible to use DA to perform finite differences on two staggered regular grids (as in the electromagnetic finite difference time domain method) ? Surrounding nodes from one grid are used to update the value in the dual grid. In addition, local manipulations need to be done on the nodal values. Thanks Rgds, Amit
SNES Convergence
Hi, I am solving a version of the Grad-Shafranov equation, F(x)=0, which but for some extra spatial dependences is similar in form to the 2D Bratu equation in snes/examples/tutorials/ex5.c. I started with the ex5.c code, introducing just enough changes to model the new system. The analytic Jacobian function appears to be correct, with a Norm of matrix ratio 1.e9(found using -snes_type test). The problem I am having is that in most cases the SNES solver does not converge to a solution x of F(x)=0. Rather, what happens is that the fnorm (obtained in a monitor function) converges to some large non-zero value, and F(x) seems to get stuck, i.e. it converges to a large non-zero result. Even though it is clearly not a solution, the -snes_converged_reason is reported as Nonlinear solve converged due to CONVERGED_TR_DELTA. The intermediate KSP steps have -ksp_converged_reason reported as Linear solve converged due to CONVERGED_STEP_LENGTH. I have been typically running with parameters -da_grid_x 100 -da_grid_y 101 -snes_converged_reason -ksp_converged_reason -snes_type tr -ksp_type gmres -snes_max_it 100. Does this sound like a familiar scenario with a familiar solution? Or can anyone point me to some documentation that describes the SNES tr and ls parameters in more detail than the manual.pdf? Or can anyone recommend the best SNES and KSP parameters for the Bratu example? Any help or advice would be greatly appreciated. Thanks, Sean Dettrick -- next part -- An HTML attachment was scrubbed... URL: http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20080107/7ec37516/attachment.htm
SNES Convergence
On Jan 7, 2008 4:35 PM, Barry Smith bsmith at mcs.anl.gov wrote: Sean, Can you run the ls version with the additional options -info - snes_monitor and send the results? Hi Barry, Attached are the gzipped results. It says something about an inconsistent rhs. Could that be the problem? And what would the rhs in question be? You could also make a run with -snes_mf_operator and see if that converges in the same way or not (this is another check that the analytic Jacobian is right/wrong). I tried that, and it seemed to do the same thing. Thanks a lot for the feedback. Sean Thanks Barry On Jan 7, 2008, at 3:30 PM, Sean Dettrick wrote: Hi, I am solving a version of the Grad-Shafranov equation, F(x)=0, which but for some extra spatial dependences is similar in form to the 2D Bratu equation in snes/examples/tutorials/ex5.c. I started with the ex5.c code, introducing just enough changes to model the new system. The analytic Jacobian function appears to be correct, with a Norm of matrix ratio 1.e9 (found using -snes_type test). The problem I am having is that in most cases the SNES solver does not converge to a solution x of F(x)=0. Rather, what happens is that the fnorm (obtained in a monitor function) converges to some large non-zero value, and F(x) seems to get stuck, i.e. it converges to a large non-zero result. Even though it is clearly not a solution, the -snes_converged_reason is reported as Nonlinear solve converged due to CONVERGED_TR_DELTA. The intermediate KSP steps have -ksp_converged_reason reported as Linear solve converged due to CONVERGED_STEP_LENGTH. I have been typically running with parameters -da_grid_x 100 -da_grid_y 101 -snes_converged_reason - ksp_converged_reason -snes_type tr -ksp_type gmres -snes_max_it 100. Does this sound like a familiar scenario with a familiar solution? Or can anyone point me to some documentation that describes the SNES tr and ls parameters in more detail than the manual.pdf? Or can anyone recommend the best SNES and KSP parameters for the Bratu example? Any help or advice would be greatly appreciated. Thanks, Sean Dettrick -- next part -- An HTML attachment was scrubbed... URL: http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20080107/e2491c32/attachment.htm -- next part -- A non-text attachment was scrubbed... Name: ls_result.gz Type: application/x-gzip Size: 2600 bytes Desc: not available URL: http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20080107/e2491c32/attachment.bin
c++ bindings, intel compiler, petsc.h conflict, petsc install?
Hi, I have a C++ MPI code that uses the C++ MPI bindings and uses PETSc as a kind of plug-in. This was compiling and running nicely before on linux with petsc-2.3.0, GCC compilers and LAM MPI. But now I can't get my code to compile on linux with petsc-2.3.3-p4, 64 bit Intel compilers, and Intel MPI. If petsc.h is included before mpi.h, then there are compile time errors: the MPI namespace is not available, due to MPICH_SKIP_MPICXX in petsc.h (not present in 2.3.0). If mpi.h is included before petsc.h, then there are link time errors: undefined referenced to Petsc functions. Am I doing something wrong, or is this a well known thing? I wonder if it could be a configuration error? Because the intel mpi installation uses non-standard directory names for the 64 bit version - bin64, include64 etc, the configuration looks a little hairy: config/configure.py \ --with-petsc-arch=intel_MPI_64_static \ --with-fortran=0 \ --with-mpi-include=/opt/intel/mpi/3.0/include64 \ --with-mpi-lib=/opt/intel/mpi/3.0/lib64/libmpi.a \ --with-mpiexec=/opt/intel/mpi-rt/3.0/bin64/mpiexec \ --with-x=no \ --with-matlab=0 \ --with-shared=0 \ --with-cc=/opt/intel/mpi/3.0/bin64/mpiicc \ --with-cxx=/opt/intel/mpi/3.0/bin64/mpiicpc \ --CXXFLAGS=-I/opt/intel/mpi/3.0/include64 \ --CFLAGS=-I/opt/intel/mpi/3.0/include64 \ --LDFLAGS=-L/opt/intel/mpi/3.0/lib64 \ --download-c-blas-lapack=1 make all test indicated that the tests were passed. Any suggestions greatly appreciated! Thanks Sean Dettrick
polar coordinate singularity
Hi, I am solving Poisson's equation Div. Grad(Phi)=rho in the (R,Theta) plane with KSP, and it is working well if I avoid the axis. But I would like to put a point at R=0. Is there a recommended way to do that? Thanks, Sean
