Re: [petsc-users] Eigenvalues always converge to zero when using slepc4py-complex

2021-10-06 Thread Lisandro Dalcin
The new builds are up. Jose added instructions in the SLEPc FAQ. Could you
try again?

Regards,

On Mon, 4 Oct 2021 at 16:13 Yelyzaveta Velizhanina 
wrote:

> I see. Thanks, much appreciated.
>
>
>
> Best regards,
>
> Yelyzaveta Velizhanina.
>
>
>
> *From: *Jose E. Roman 
> *Date: *Monday, 4 October 2021 at 15:10
> *To: *Yelyzaveta Velizhanina 
> *Cc: *petsc-users@mcs.anl.gov 
> *Subject: *Re: [petsc-users] Eigenvalues always converge to zero when
> using slepc4py-complex
>
> Conda supports complex scalars for petsc4py. However, this is not
> implemented in slepc4py. Lisandro is trying to get this fixed, so if no
> issues arise this will be available in a couple of days, with
> slepc4py-3.16.0.
>
> Jose
>
>
> > El 3 oct 2021, a las 22:45, Yelyzaveta Velizhanina <
> velizhani...@gmail.com> escribió:
> >
> > Dear all,
> >
> > I am having a problem to get EPS run properly with PETSc and SLEPc build
> with scalar_value=complex. I am using petsc4py and slepc4py. Installed
> everything, including PETSc and SLEPc, with conda. While real scalar value
> build works well, when using the complex one, all the eigenvalues always
> converge to 0 for any matrix and any solver. I’ve tried running examples
> given in this repo https://github.com/myousefi2016/slepc4py as well -
> same outcome, only zero eigenvalues. I am running MacOSX BigSur.
> >
> > Will appreciate any help,
> >
> > Best regards,
> > Yelyzaveta Velizhanina.
>
-- 
Lisandro Dalcin

Senior Research Scientist
Extreme Computing Research Center (ECRC)
King Abdullah University of Science and Technology (KAUST)
http://ecrc.kaust.edu.sa/


Re: [petsc-users] headsup: switch git default branch from 'master' to 'main'

2021-02-23 Thread Lisandro Dalcin
May the force be with you, Satish.

On Tue, 23 Feb 2021 at 20:19, Satish Balay via petsc-users <
petsc-users@mcs.anl.gov> wrote:

> All,
>
> This is a heads-up, we are to switch the default branch in petsc git
> repo from 'master' to 'main'
>
> [Will plan to do the switch on friday the 26th]
>
> We've previously switched 'maint' branch to 'release' before 3.14
> release - and this change (to 'main') is the next step in this direction.
>
> Satish
>
>

-- 
Lisandro Dalcin

Senior Research Scientist
Extreme Computing Research Center (ECRC)
King Abdullah University of Science and Technology (KAUST)
http://ecrc.kaust.edu.sa/


Re: [petsc-users] Python version needed for internal scripts

2020-10-14 Thread Lisandro Dalcin
On Wed, 14 Oct 2020 at 14:12, Lisandro Dalcin  wrote:

> Disclaimer: I'm not an Ubuntu user. Google and inform yourself.
>

Just in case... What I'm trying to say is that I do not know Ubuntu very
well, much less its  Py2 -> Py3 transition details, and then you should try
to find a more authoritative source than my occasional comments about the
proper way to do things.


-- 
Lisandro Dalcin

Senior Research Scientist
Extreme Computing Research Center (ECRC)
King Abdullah University of Science and Technology (KAUST)
http://ecrc.kaust.edu.sa/


Re: [petsc-users] Python version needed for internal scripts

2020-10-14 Thread Lisandro Dalcin
On Wed, 14 Oct 2020 at 13:01, Nicolas Barral <
nicolas.bar...@math.u-bordeaux.fr> wrote:

>
> 'python' is usually an alias for python2, so making it point at python3
> seems a bit dangerous. Yet, python2 was removed from recent Ubuntus and
> maybe others, and if I have no python2 installed, and no 'python' alias,
> I have to manually edit all the scripts.
>
>
apt install python-is-python3

and you should get the alias python -> python3 in /usr/bin

Disclaimer: I'm not an Ubuntu user. Google and inform yourself.

-- 
Lisandro Dalcin

Senior Research Scientist
Extreme Computing Research Center (ECRC)
King Abdullah University of Science and Technology (KAUST)
http://ecrc.kaust.edu.sa/


Re: [petsc-users] compiling PETSc using c++ compiler

2020-10-07 Thread Lisandro Dalcin
On Wed, 7 Oct 2020 at 02:45, Sam Guo  wrote:

>I want to experiment c++ compiler to see if I can compile real and
> complex versions to different symbols (I've created another thread for this
> topic.) but it seems c++ compiler does not help and I still get same
> symbols.
>
>
Indeed, that is not going to make it. Or you have to change the definition
of PETSC_EXTERN, such that it does not use `extern "C"`.  Even with that,
you may have trouble with EXTERN_C_BEGIN/END macros.

-- 
Lisandro Dalcin

Senior Research Scientist
Extreme Computing Research Center (ECRC)
King Abdullah University of Science and Technology (KAUST)
http://ecrc.kaust.edu.sa/


Re: [petsc-users] HDF5 and ParaView

2020-07-13 Thread Lisandro Dalcin
On Mon, 13 Jul 2020 at 14:14, Mark Adams  wrote:

>
>>
>> I would probably move out of VTK files in favor of something else if I
>> had a way to encode VTK's (the library, not the file format) high-order
>> Lagrange elements.
>> Actually, I'm toying with dumping files with PETSc's raw binary I/O with
>> MPI, and writing a proper ParaView plugin in Python to read the data.
>>
>
> I'd love a high order viewer too.
>

Well, ParaView can actually do it. You need to open a vtu file with
high-order Lagrange elements, and then look for "Nonlinear subdivision
level" in the GUI and bump it up. The main issue is that the algorithm is
not "adaptive" and you may need to bump too much, then the thing becomes
quite slow.

Another one that does things very well is GLVis, although not with all the
features of the ParaView/VisIt+VTK combo. Stefano put a lot of effort
integrating GLVIs into PETSc, but maybe the integration with PetscFE and
DMPlex may require some extra work.

-- 
Lisandro Dalcin

Research Scientist
Extreme Computing Research Center (ECRC)
King Abdullah University of Science and Technology (KAUST)
http://ecrc.kaust.edu.sa/


Re: [petsc-users] HDF5 and ParaView

2020-07-13 Thread Lisandro Dalcin
On Fri, 10 Jul 2020 at 15:37, Matthew Knepley  wrote:

>
> That is true. Do they also get rid of the single mesh and single timstep
> requirements? HDF5+XDMF makes it much easier since
> we can put multiple meshes and timesteps in one file.
>

Well, I use ParaView's *.pvd files for that, and I dump each timestep to
its own *.vtu file in a folder to pack the files. But you still have a
point, the format is indeed restricting.
Why do you consider it so important to put multiple meshes and timesteps in
one file? It is just that you hate to have so many files scattered around?
Or something deeper?
I do hate the fact that the VTK formats (either legacy or XML) do not allow
you to dump a single mesh to be reused for multiple timestep.

I would probably move out of VTK files in favor of something else if I had
a way to encode VTK's (the library, not the file format) high-order
Lagrange elements.
Actually, I'm toying with dumping files with PETSc's raw binary I/O with
MPI, and writing a proper ParaView plugin in Python to read the data.


-- 
Lisandro Dalcin

Research Scientist
Extreme Computing Research Center (ECRC)
King Abdullah University of Science and Technology (KAUST)
http://ecrc.kaust.edu.sa/


Re: [petsc-users] HDF5 and ParaView

2020-07-10 Thread Lisandro Dalcin
On Thu, 9 Jul 2020 at 23:11, Matthew Knepley  wrote:

> On Thu, Jul 9, 2020 at 3:28 PM David Scott  wrote:
>
>> Matt,
>>
>> I'll have a go at writing Xdmf files.
>>
>> Am I right in thinking that VTK files are written sequentially?
>>
>
> Yes, it is a rightly despised format, suitable only for beggars and serial
> jobs :)
>
>
By "rightly despised format", I hope you are talking about the legacy VTK
files (the traditional .vtk extension).
The new XML formats with  are quite easy to
write in parallel with MPI/IO (at least *.vtu files, which is what I'm
using).

-- 
Lisandro Dalcin

Research Scientist
Extreme Computing Research Center (ECRC)
King Abdullah University of Science and Technology (KAUST)
http://ecrc.kaust.edu.sa/


Re: [petsc-users] Terminating a process running petsc via petsc4py without mpi_abort

2020-06-04 Thread Lisandro Dalcin
(1) You can use PETSc.Sys.pushErrorHandler("abort"), but it will not help
you. What you really need is to override PETSc's default signal handling

(2) While it is true that PETSc overrides the signal handler, you can
override it again from python after from petsc4py import PETSc.

For implementing (2), maybe you should try sending SIGINT and not SIGTERM,
such that you can do the following.

from petsc4py import PETSc

import signal
signal.signal(signal.SIGINT, signal.default_int_handler)

...

if __name__ == "__main__":
try:
main()
except KeyboardInterrupt: # Triggered if Ctrl+C or signaled with SIGINT
... # do cleanup if needed

Otherwise, you just need  signal.signal(signal.SIGINT, signal.SIG_DFL)


PS: I'm not in favor of changing current PETSc's signal handling behavior.
This particular issue is fixable with two lines of Python code:

from signal import signal, SIGINT, SIG_DFL
signal(SIGINT, SIG_DFL)



On Thu, 4 Jun 2020 at 23:39, Hudson, Stephen Tobias P 
wrote:

> Lisandro,
>
> I don't see an interface to set this through petsc4py. Is it possible?
>
> Thanks,
> Steve
> --
> *From:* Hudson, Stephen Tobias P 
> *Sent:* Thursday, June 4, 2020 2:47 PM
> *To:* Balay, Satish 
> *Cc:* petsc-users@mcs.anl.gov ; Lisandro Dalcin <
> dalc...@gmail.com>
> *Subject:* Re: Terminating a process running petsc via petsc4py without
> mpi_abort
>
> Sounds good. I will have a look at how to set this through petsc4py.
>
> Thanks
> Steve
> --
> *From:* Satish Balay 
> *Sent:* Thursday, June 4, 2020 2:32 PM
> *To:* Hudson, Stephen Tobias P 
> *Cc:* petsc-users@mcs.anl.gov ; Lisandro Dalcin <
> dalc...@gmail.com>
> *Subject:* Re: Terminating a process running petsc via petsc4py without
> mpi_abort
>
> I don't completely understand the issue here. How is sequential run
> different than parallel run?
>
> In both cases - a PetscErrorHandler is likely getting invoked. One can
> change this behavior with:
>
>
> https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Sys/PetscPushErrorHandler.html
>
> And there are a few default error handlers to choose
>
>
> PETSC_EXTERN PetscErrorCode PetscTraceBackErrorHandler(MPI_Comm,int,const
> char*,const char*,PetscErrorCode,PetscErrorType,const char*,void*);
> PETSC_EXTERN PetscErrorCode PetscIgnoreErrorHandler(MPI_Comm,int,const
> char*,const char*,PetscErrorCode,PetscErrorType,const char*,void*);
> PETSC_EXTERN PetscErrorCode
> PetscEmacsClientErrorHandler(MPI_Comm,int,const char*,const
> char*,PetscErrorCode,PetscErrorType,const char*,void*);
> PETSC_EXTERN PetscErrorCode PetscMPIAbortErrorHandler(MPI_Comm,int,const
> char*,const char*,PetscErrorCode,PetscErrorType,const char*,void*);
> PETSC_EXTERN PetscErrorCode PetscAbortErrorHandler(MPI_Comm,int,const
> char*,const char*,PetscErrorCode,PetscErrorType,const char*,void*);
> PETSC_EXTERN PetscErrorCode
> PetscAttachDebuggerErrorHandler(MPI_Comm,int,const char*,const
> char*,PetscErrorCode,PetscErrorType,const char*,void*);
> PETSC_EXTERN PetscErrorCode PetscReturnErrorHandler(MPI_Comm,int,const
> char*,const char*,PetscErrorCode,PetscErrorType,const char*,void*);
>
> Some of the are accessible via command line option. for ex:
> -on_error_abort or -on_error_mpiabort
>
> Or perhaps you want to completely disable error handler with:
> -no_signal_handler
>
> cc: petsc-users
>
> Satish
>
> On Thu, 4 Jun 2020, Hudson, Stephen Tobias P wrote:
>
> > Satish,
> >
> > We are having issues caused by MPI_abort getting called when we try to
> terminate a sub-process running petsc4py. Ideally we would always use a
> serial build of petsc/petsc4py in this mode, but many users will have a
> parallel build. We need to be able to send a terminate signal that just
> kills the process.
> >
> > Is there a way to turn off the mpi_abort?
> >
> > Thanks,
> >
> > Steve
> >
> >
>
>

-- 
Lisandro Dalcin

Research Scientist
Extreme Computing Research Center (ECRC)
King Abdullah University of Science and Technology (KAUST)
http://ecrc.kaust.edu.sa/


Re: [petsc-users] Resuming analysis by importing Trajectory Data from external source

2020-05-10 Thread Lisandro Dalcin
On Sat, 9 May 2020 at 23:44, Mohammed Ashour  wrote:

> I'm using TSALPHA1, as in ierr = TSSetType(ts,TSALPHA1); CHKERRQ(ierr);
>
>
Well, TSALPHA1 uses some vectors keeping intenal state, as you know,
generalized-\alpha is not really a one-step method. That intenal state
should be saved and loaded if you want to implement a truly correct
checkpoint/restart. But PETSc does not currently provide access to this
vector in its public API, so you will have to hack things around. IIRC,
just saving and restoring the internal V1 vector is all what you need,
maybe also set ts->steprestart to FALSE and properly set the ts->steps
counter.

Or you can just ignore these detail, and restart the easy way from just the
saved solution, at the price of a very small "glitch" at the restart time.
But note that this glitch just means that you temporarily change the scheme
to reinitialize, and it is exactly the same thing the implementation does
at the initial time to initialize the method (as you surely know, the
method is also not trully self-starting). The restarting procesure as
implemented in TSALPHA1 is not in the literature, it is of my own cooking
but based on rather trivial relations performing two steps with dt/2 time
step size, the fist of those inner steps using backward-Euler.


-- 
Lisandro Dalcin

Research Scientist
Extreme Computing Research Center (ECRC)
King Abdullah University of Science and Technology (KAUST)
http://ecrc.kaust.edu.sa/


Re: [petsc-users] Resuming analysis by importing Trajectory Data from external source

2020-05-09 Thread Lisandro Dalcin
On Sat, 9 May 2020 at 22:48, Mohammed Ashour  wrote:

> Oh, thanks for the tip.
> I can do that. Aside from the solution vector, step number, and time, are
> there any hidden details in TS I need to dump that would be essential for
> the next run?
>
>
This really depends on the timestepper.  What TS type are you using?

-- 
Lisandro Dalcin

Research Scientist
Extreme Computing Research Center (ECRC)
King Abdullah University of Science and Technology (KAUST)
http://ecrc.kaust.edu.sa/


Re: [petsc-users] [EXTERNAL] Re: Alternate mesh partitioners to METIS/ParMETIS

2020-05-07 Thread Lisandro Dalcin
On Thu, 7 May 2020 at 17:53, Hammond, Glenn E via petsc-users <
petsc-users@mcs.anl.gov> wrote:

> Jed,
>
> We call MatMeshToCellGraph() to generate the dual matrix.  This function
> relies upon ParMETIS.  Can you point me to a similar function that does not
> require ParMETIS ( e.g. for PTScotch)?
>
>
Maybe a starting point:

DMPlexCreateFromCellList(comm, ..., );
PetscSectionCreate(comm, );
DMPlexGetPartitioner(dm, );
PetscPartitionerSetType(partitioner, PETSCPARTITIONERPTSCOTCH);
PetscPartitionerDMPlexPartition(partitioner, dm, NULL, cellPartSection,
);

One minor annoyance is that the first call will need the vertex
coordinates.  Matthew, any better way?


-- 
Lisandro Dalcin

Research Scientist
Extreme Computing Research Center (ECRC)
King Abdullah University of Science and Technology (KAUST)
http://ecrc.kaust.edu.sa/


Re: [petsc-users] PetscObjectCompose error

2020-04-11 Thread Lisandro Dalcin
On Sun, 12 Apr 2020 at 01:50, Mark Adams  wrote:

> I am trying to compose a ISColoring to a Mat. THe code works, I know JacP
> and iscoloring are valid Mat and ISColoring. I have this:
>
> ierr =
> ((PetscObject)JacP,"coloring",(PetscObject)iscoloring);CHKERRQ(ierr);
>
> But it says my ISColoring is not valid. Any suggestions?
>
>
ISColoring is not a PetscObject, you cannot compose it.

struct _n_ISColoring {
  PetscIntrefct;
  PetscIntn;/* number of colors */
  IS  *is;  /* for each color indicates columns */
  MPI_Commcomm;
  ISColoringValue *colors;  /* for each column indicates color */
  PetscIntN;/* number of columns */
  ISColoringType  ctype;
  PetscBool   allocated;
};

I guess your best option for now is to store it in a
PetscContainer object with a destroy callback, and then compose the
container.



-- 
Lisandro Dalcin

Research Scientist
Extreme Computing Research Center (ECRC)
King Abdullah University of Science and Technology (KAUST)
http://ecrc.kaust.edu.sa/


Re: [petsc-users] How to set an initial guess for TS

2020-04-03 Thread Lisandro Dalcin
On Fri, 3 Apr 2020 at 18:26, Fande Kong  wrote:

> Hi All,
>
> TSSetSolution will set an initial condition for the current TSSolve().
> What should I do if I want to set an initial guess for the current solution
> that is different from the initial condition?  The initial guess is
> supposed to be really close to the current solution, and then will
> accelerate my solver.
>
> In other words, TSSetSolution will set "U_{n-1}", and now we call TSSolve
> to figure out "U_{n}". If I know something about "U_{n}", and I want to set
> "\bar{U}_{n}" as the initial guess of "U_{n}" when computing "U_{n}".
>
>
IMHO, the only reliable solution would be to implement
SNESSet{Pre|Post}Solve(snes, {Pre|Post}Solve, ctx) to set callback routines
{Pre|Post}Solve(snes,b,x,ctx) allowing users to modify the solution vector
'x' in place (but no 'b', of course).
We already have that feature in KSP, why not SNES which is an inner loop in
TS ?


-- 
Lisandro Dalcin

Research Scientist
Extreme Computing Research Center (ECRC)
King Abdullah University of Science and Technology (KAUST)
http://ecrc.kaust.edu.sa/


Re: [petsc-users] Error Installing Petsc4Py

2020-04-02 Thread Lisandro Dalcin
Please note that I have not released yet petsc4py-3.13, so pip install will
not work with latest PETSc release 3.13

On Thu, 2 Apr 2020 at 20:48, Satish Balay via petsc-users <
petsc-users@mcs.anl.gov> wrote:

> Which version of petsc and which version of petsc4py is this?
>
> Lisandro can help with pip issues.
>
> Just want to mention an alternate mode that you might want to try:
>
> --download-petsc4py option to petsc configure
>
> Satish
>
> On Thu, 2 Apr 2020, Kyle Michael Booker wrote:
>
> > Hello,
> >
> >
> > I've been trying to install petsc4py for the past several days and have
> had the following errors installing:
> >
> >
> > "Command "/usr/bin/python -u -c "import setuptools,
> tokenize;__file__='/tmp/pip-build-0eWEiE/petsc4py/setup.py';f=getattr(tokenize,
> 'open', open)(__file__);code=f.read().replace('\r\n',
> '\n');f.close();exec(compile(code, __file__, 'exec'))" install --record
> /tmp/pip-oUZMRL-record/install-record.txt
> --single-version-externally-managed --compile --user --prefix=" failed with
> error code 1 in /tmp/pip-build-0eWEiE/petsc4py/"
> >
> >
> > I currently have petsc installed correctly as well as mpich, it is just
> getting petsc4py installed that I have problems with. Any help would be
> greatly appreciated.
> >
> >
> > Sincerely,
> >
> >
> > Kyle
> >
>
>

-- 
Lisandro Dalcin

Research Scientist
Extreme Computing Research Center (ECRC)
King Abdullah University of Science and Technology (KAUST)
http://ecrc.kaust.edu.sa/


Re: [petsc-users] [petsc4py] Assembly fails

2020-03-27 Thread Lisandro Dalcin
On Fri, 27 Mar 2020 at 17:10, Matthew Knepley  wrote:

> On Fri, Mar 27, 2020 at 3:31 AM Alejandro Aragon - 3ME <
> a.m.ara...@tudelft.nl> wrote:
>
>> Dear Matthew,
>>
>> Thanks for your email. I have attached the python code that reproduces
>> the following error in my computer:
>>
>
> I think I see the problem. There were changes in DM in order to support
> fields which only occupy part of the domain.
> Now you need to tell the DM about the fields before it builds a Section. I
> think in your code, you only need
>
>   f = PetscContainer()
>   f.setName("potential")
>   dm.addField(field = f)
>
>
Except that petsc4py do not expose a Container class :-(

@Alejandro, PETSc 3.13 is about to be released. Once that is done, we have
a small windows to make a few enhancements in petsc4py before releasing the
matching petsc4py-3.13 that you can use to upgrade your code. Would that
work for you?

-- 
Lisandro Dalcin

Research Scientist
Extreme Computing Research Center (ECRC)
King Abdullah University of Science and Technology (KAUST)
http://ecrc.kaust.edu.sa/


Re: [petsc-users] [petsc4py] Assembly fails

2020-03-26 Thread Lisandro Dalcin
On Wed, 25 Mar 2020 at 19:38, Matthew Knepley  wrote:

>
> It looks like you have an inconsistent build, or a memory overwrite. Since
> you are in Python, I suspect the former. Can you build
> PETSc from scratch and try this? Does it work in serial? Can you send a
> small code that reproduces this?
>
>
Well, I think Alejandro installed PETSc/petsc4py via pip as per my
instructions. Can you confirm, Alejandro?



-- 
Lisandro Dalcin

Research Scientist
Extreme Computing Research Center (ECRC)
King Abdullah University of Science and Technology (KAUST)
http://ecrc.kaust.edu.sa/


Re: [petsc-users] KSPComputeOperator in petsc4py

2020-01-05 Thread Lisandro Dalcin
On Thu, 26 Dec 2019 at 21:36, Mark Cunningham <
mark.cunning...@ariacoustics.com> wrote:

> I have an application written in python and using petsc4py.  I would like
> to study the effects of different preconditioners on my eigenvalue spectrum
> and thought to use KSPComputeOperator to obtain the preconditioned matrix.
> This, apparently, is not disclosed through the petsc4py implementation.
>

Could you share the modifications you implemented? This way I could quickly
add them to petsc4py to make them available in next release.


> I found KSP.pyx, where I believe that the KSP object is defined, and added
> a definition for the function but, after a pip install,
>

If you "cd" to the top level petsc4py source directory with your
modifications, then you may need

$ pip install --no-cache-dir . # note the final dot "."

otherwise pip may just reinstall petsc4py from a previously built wheel
file stored in the pip cache.

-- 
Lisandro Dalcin

Research Scientist
Extreme Computing Research Center (ECRC)
King Abdullah University of Science and Technology (KAUST)
http://ecrc.kaust.edu.sa/


Re: [petsc-users] RAW binary write

2018-12-05 Thread Lisandro Dalcin via petsc-users
On Tue, 4 Dec 2018 at 19:53, Matthew Knepley  wrote:

> On Tue, Dec 4, 2018 at 11:11 AM Lisandro Dalcin via petsc-users <
> petsc-users@mcs.anl.gov> wrote:
>
>> PetscByteSwap ? But be careful, you should (1) vec get array, byteswap,
>> vec restore array, (2) view, (3) vec get array, byteswap, vec restore
>> array. Of course, if you do not care about the final Vec contents, you can
>> just omit step (3) and let it be, let it be...
>>
>> PS: Eventually, PETSc binary viewers should gain an API like
>> PetscViewerBinarySetByteOrder(viewer, PETSC_LITTLE_ENDIAN) ...
>>
>
> Note that this enables exactly what we wanted to avoid, namely binary
> files that an incompatible with some installations.
>
>
1) As a confessed libertarian (not to be confused with liberal), IMHO, we
should not prevent power users to do something if they explicitly ask for
it. Maybe someone has to dump a binary file in little endian with no
headers, so that some other crappy software  that do not care about
endianness can happily read the file and do something with it. Why would we
not add the feature (at least if someone contributes it)?

2) The whole native binary I/O in PETSc is messy.

3) In an ideal word, every PETSc binary file should start with a header
8bitEndianFlag,sizeof(PetscInt),sizeof(PetscReal),sizeof(PetscScalar), And
maybe readers should handle gracefully binary files produced by PETSc
builds with different integer/real sizes.

-- 
Lisandro Dalcin

Research Scientist
Computer, Electrical and Mathematical Sciences & Engineering (CEMSE)
Extreme Computing Research Center (ECRC)
King Abdullah University of Science and Technology (KAUST)
http://ecrc.kaust.edu.sa/

4700 King Abdullah University of Science and Technology
al-Khawarizmi Bldg (Bldg 1), Office # 0109
Thuwal 23955-6900, Kingdom of Saudi Arabia
http://www.kaust.edu.sa

Office Phone: +966 12 808-0459


Re: [petsc-users] RAW binary write

2018-12-02 Thread Lisandro Dalcin via petsc-users
Use PetscViewerBinarySetSkipHeader(viewer, PETSC_TRUE) before VecView().

On Thu, 29 Nov 2018 at 18:50, Sal Am via petsc-users <
petsc-users@mcs.anl.gov> wrote:

> Is there a way to write the solution from the system Ax=b in raw binary
> instead of PETSc binary format?
>
> Currently I am doing:
>   ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,
> "../../python/petscpy/Vector_x_petsc.dat", FILE_MODE_WRITE,
> );CHKERRQ(ierr);
>   ierr = VecView(x,viewer);CHKERRQ(ierr);CHKERRQ(ierr);
>   ierr = PetscViewerDestroy();CHKERRQ(ierr);
>
> And then use PetscBinaryIO to read it back and save it using
> write('newformat', 'wb') to get to raw... however this approach is not good
> it seems as there are some troubles with little/big endian when using the
> resulting converted file on other systems for post-processing.
>
> Thanks,
>


-- 
Lisandro Dalcin

Research Scientist
Computer, Electrical and Mathematical Sciences & Engineering (CEMSE)
Extreme Computing Research Center (ECRC)
King Abdullah University of Science and Technology (KAUST)
http://ecrc.kaust.edu.sa/

4700 King Abdullah University of Science and Technology
al-Khawarizmi Bldg (Bldg 1), Office # 0109
Thuwal 23955-6900, Kingdom of Saudi Arabia
http://www.kaust.edu.sa

Office Phone: +966 12 808-0459


Re: [petsc-users] Extending petsc4py

2018-07-11 Thread Lisandro Dalcin
On Tue, 10 Jul 2018 at 20:17, Ingo Gaertner 
wrote:

> Hello,
> can I find somewhere an example how to wrap my own PETSc-based code so
> that it can be used together with petsc4py?
>
> Let's assume I have a C function, which I want to access from python:
> Vec createMyVec();
>
> I have written a SWIG interface file to generate a python wrapper for this
> function. The wrapper works, but I cannot use the return value as input for
> mpi4py routines that expect a Vec, because mpi4py's Vec type is different
> from the SWIG generated Vec type in my python wrapper. Any idea, how I can
> create mpi4py-compatible types?
>
> If the mpi4py wrapping technique is described somewhere, this may help as
> well, although some working example would be perfect.
>
> As a final remark, I'd like to keep my project decoupled from mpi4py.
> Therefore, including my extensions in the build process of mpi4py is not
> quite what I am looking for, although I could live with this as a
> workaround if no other solution is possible.
> (I am using petsc-3.9.1 built from source and mpi4py installed using pip
> on Ubuntu 18.04.)
>
>
Sorry, but I'm quite confused. Why are you talking about mpi4py? There is
no "Vec" in mpi4py. The only type intersection between petsc4py and mpi4py
is "Comm" for MPI communicators, in that case you should use the following
to make the conversions:

mpi4py_comm = petsc4py_comm.tompi4py() # mpi4py <-- petsc4py
petsc4py_comm = PETSc.Comm(mpi4py_comm) # petsc4py <-- mpi4py

About SWIG, are you using typemaps written by yourself, or the typemaps
provided by petsc4py? Have you looked at "demo/wrap-swig" in petsc4py
sources?

-- 
Lisandro Dalcin

Research Scientist
Computer, Electrical and Mathematical Sciences & Engineering (CEMSE)
Extreme Computing Research Center (ECRC)
King Abdullah University of Science and Technology (KAUST)
http://ecrc.kaust.edu.sa/

4700 King Abdullah University of Science and Technology
al-Khawarizmi Bldg (Bldg 1), Office # 0109
Thuwal 23955-6900, Kingdom of Saudi Arabia
http://www.kaust.edu.sa

Office Phone: +966 12 808-0459


Re: [petsc-users] "snes_type test" is gone?

2018-06-30 Thread Lisandro Dalcin
./example.exe -snes_test_jacobian -snes_type ksponly -ksp_type preonly
-pc_type none -snes_convergence_test skip
On Sat, 30 Jun 2018 at 01:06, Alexander Lindsay
 wrote:
>
> Is there a way to mimic the behavior of `-snes_type test` with new PETSc? 
> E.g. don't attempt to perform any solves? They're stupid but we have a bunch 
> of tests in MOOSE that are set-up only to test the Jacobian, and any attempts 
> to actually solve the system are disastrous. I can hack around this by doing 
> things like setting the convergence reason to `SNES_CONVERGED_ITS` and 
> similarly for ksp, making sure I use `-pc_type gamg` since there are missing 
> Jacobian diagonal entries, etc. but this doesn't seem ideal.
>
> Alex
>
> On Mon, Jun 4, 2018 at 4:26 PM, Kong, Fande  wrote:
>>
>> Thanks, Hong,
>>
>> I see. It is better if "-snes_type test" did not exist at the first place.
>>
>>
>> Fande,
>>
>> On Mon, Jun 4, 2018 at 4:01 PM, Zhang, Hong  wrote:
>>>
>>>
>>>
>>> > On Jun 4, 2018, at 4:59 PM, Zhang, Hong  wrote:
>>> >
>>> > -snes_type has been removed. We can just use -snes_test_jacobian instead. 
>>> > Note that the test is done every time the Jacobian is computed.
>>>
>>> It was meant to be "-snes_type test".
>>>
>>> > Hong (Mr.)
>>> >
>>> >> On Jun 4, 2018, at 3:27 PM, Kong, Fande  wrote:
>>> >>
>>> >> Hi PETSc Team,
>>> >>
>>> >> I was wondering if "snes_type test" has been gone? Quite a few MOOSE 
>>> >> users use this option to test their Jacobian matrices.
>>> >>
>>> >> If it is gone, any reason?
>>> >>
>>> >> Fande,
>>> >
>>>
>>
>


-- 
Lisandro Dalcin

Research Scientist
Computer, Electrical and Mathematical Sciences & Engineering (CEMSE)
Extreme Computing Research Center (ECRC)
King Abdullah University of Science and Technology (KAUST)
http://ecrc.kaust.edu.sa/

4700 King Abdullah University of Science and Technology
al-Khawarizmi Bldg (Bldg 1), Office # 0109
Thuwal 23955-6900, Kingdom of Saudi Arabia
http://www.kaust.edu.sa

Office Phone: +966 12 808-0459


Re: [petsc-users] MATIS and DMDA in petsc4py

2018-04-26 Thread Lisandro Dalcin
On 26 April 2018 at 11:05, Loic Gouarin <loic.goua...@math.u-psud.fr> wrote:
> The problem is solved with the new 3.9.0 release on conda-forge.
>
> Thanks Lisandro !!
>

Well, you are most welcome, but I don't remember doing anything in
particular to get that fixed :-)

-- 
Lisandro Dalcin

Research Scientist
Computer, Electrical and Mathematical Sciences & Engineering (CEMSE)
Extreme Computing Research Center (ECRC)
King Abdullah University of Science and Technology (KAUST)
http://ecrc.kaust.edu.sa/

4700 King Abdullah University of Science and Technology
al-Khawarizmi Bldg (Bldg 1), Office # 0109
Thuwal 23955-6900, Kingdom of Saudi Arabia
http://www.kaust.edu.sa

Office Phone: +966 12 808-0459


Re: [petsc-users] petsc4py: reuse setup for multiple solver calls?

2018-04-06 Thread Lisandro Dalcin
On 6 April 2018 at 08:48, Robert Speck <r.sp...@fz-juelich.de> wrote:
> Thank you for your answer! Please see below for comments/questions.
>
> On 05.04.18 12:53, Matthew Knepley wrote:
>> On Thu, Apr 5, 2018 at 6:39 AM, Robert Speck <r.sp...@fz-juelich.de

>>
>> For nonlinear solves which stay the same size, you do nothing.
>
> "nothing" in terms of "nothing you can do" or "nothing you have to do"?
>

I would say "nothing you have to do". Nonlinear solvers work by just
setting two callback routines to compute the residual and the
Jacobian. The user-defined residual routine receives a vector you have
to put values in, same for the Jacobian routine, you have to put
values in the (already preallocated) matrix that is reused on each
newton step, and on each successive solve. SNES will take care of
everything else (in particular, updating the KSP linear solver, and
refreshing the preconditioner, always trying to reuse
already-preallocated internal data structures as much as possible).
All what what I said for SNES applies to TS, in case you ever want to
give it a chance.



-- 
Lisandro Dalcin

Research Scientist
Computer, Electrical and Mathematical Sciences & Engineering (CEMSE)
Extreme Computing Research Center (ECRC)
King Abdullah University of Science and Technology (KAUST)
http://ecrc.kaust.edu.sa/

4700 King Abdullah University of Science and Technology
al-Khawarizmi Bldg (Bldg 1), Office # 0109
Thuwal 23955-6900, Kingdom of Saudi Arabia
http://www.kaust.edu.sa

Office Phone: +966 12 808-0459


Re: [petsc-users] Transform scipy sparse to partioned, parallel petsc matrix in PETSc4py

2018-02-13 Thread Lisandro Dalcin
On 11 February 2018 at 20:35, Jan Grießer <griesser@googlemail.com> wrote:
> Hey,
>
> i have a precomputed scipy sparse matrix for which I want to solve the
> eigenvalue problem for a matrix of size 35000x35000. I don´t really get how
> to parallelize this problem correctly.
> Similar to another
> thread(https://lists.mcs.anl.gov/mailman/htdig/petsc-users/2013-August/018501.html)
> I tried the following code:
>
>
>
> D = scipy.sparse.load_npz("sparse.npz")
>
> B = D.tocsr()
>
>
>
> # Construct the matrix Ds in parallel
>
> Ds = PETSc.Mat().create()
>
> Ds.setSizes(CSRmatrix.shape)
>
> Ds.assemble()
>

Use DS.setUp()

>
>
> # Fill the matrix
>
> rstart, rend = Ds.getOwnershipRange()
>
> csr = (
>
> B.indptr[rstart:rend+1] - B.indptr[rstart],
>
> B.indices[B.indptr[rstart]:B.indptr[rend]],
>
> B.data[B.indptr[rstart]:B.indptr[rend]]
>
> )
>

This looks just fine

>
>
> Ds = PETSc.Mat().createAIJ(size=CSRmatrix.shape, csr=csr)
>
> Ds.assemble()
>

I think you don't need to assemble here.

>
>
> # Solve the eigenvalue problem
>
> solve_eigensystem(Ds)
>
>
>
> This code works for 1 processor with mpiexec –n 1 python example.py, however
> for increasing number of processors it appears as if al processors try to
> solve the overall problem instead of splitting it into blocks and solve for
> a subset of eigenvalues and eigenvectors.
> Why is this the case or did I miss something?
>

I guess you are using `mpiexec` from a different MPI implementation
than the one you used to build PETSc and petsc4py.

-- 
Lisandro Dalcin

Research Scientist
Computer, Electrical and Mathematical Sciences & Engineering (CEMSE)
Extreme Computing Research Center (ECRC)
King Abdullah University of Science and Technology (KAUST)
http://ecrc.kaust.edu.sa/

4700 King Abdullah University of Science and Technology
al-Khawarizmi Bldg (Bldg 1), Office # 0109
Thuwal 23955-6900, Kingdom of Saudi Arabia
http://www.kaust.edu.sa

Office Phone: +966 12 808-0459


Re: [petsc-users] petsc4py f2py

2018-01-18 Thread Lisandro Dalcin
Line 71 is likely wrong, too, you likely meant:

do k=grd%zs,grd%ze

On 18 January 2018 at 21:59, Lisandro Dalcin <dalc...@gmail.com> wrote:
> On 18 January 2018 at 18:56, Aurelien Ponte <aurelien.po...@ifremer.fr> wrote:
>> The code is found here:
>> https://github.com/apatlpo/qgsolver/tree/fwrapping/qgsolver
>
> Your code is not valid Fortran, you cannot use brackets [] for
> indexing in lines 80-90, you have to use parentheses ().
>
> --
> Lisandro Dalcin
> 
> Research Scientist
> Computer, Electrical and Mathematical Sciences & Engineering (CEMSE)
> Extreme Computing Research Center (ECRC)
> King Abdullah University of Science and Technology (KAUST)
> http://ecrc.kaust.edu.sa/
>
> 4700 King Abdullah University of Science and Technology
> al-Khawarizmi Bldg (Bldg 1), Office # 0109
> Thuwal 23955-6900, Kingdom of Saudi Arabia
> http://www.kaust.edu.sa
>
> Office Phone: +966 12 808-0459



-- 
Lisandro Dalcin

Research Scientist
Computer, Electrical and Mathematical Sciences & Engineering (CEMSE)
Extreme Computing Research Center (ECRC)
King Abdullah University of Science and Technology (KAUST)
http://ecrc.kaust.edu.sa/

4700 King Abdullah University of Science and Technology
al-Khawarizmi Bldg (Bldg 1), Office # 0109
Thuwal 23955-6900, Kingdom of Saudi Arabia
http://www.kaust.edu.sa

Office Phone: +966 12 808-0459


Re: [petsc-users] petsc4py f2py

2018-01-18 Thread Lisandro Dalcin
On 18 January 2018 at 18:56, Aurelien Ponte <aurelien.po...@ifremer.fr> wrote:
> The code is found here:
> https://github.com/apatlpo/qgsolver/tree/fwrapping/qgsolver

Your code is not valid Fortran, you cannot use brackets [] for
indexing in lines 80-90, you have to use parentheses ().

-- 
Lisandro Dalcin

Research Scientist
Computer, Electrical and Mathematical Sciences & Engineering (CEMSE)
Extreme Computing Research Center (ECRC)
King Abdullah University of Science and Technology (KAUST)
http://ecrc.kaust.edu.sa/

4700 King Abdullah University of Science and Technology
al-Khawarizmi Bldg (Bldg 1), Office # 0109
Thuwal 23955-6900, Kingdom of Saudi Arabia
http://www.kaust.edu.sa

Office Phone: +966 12 808-0459


Re: [petsc-users] petsc4py, cython

2018-01-18 Thread Lisandro Dalcin
In Cython code using C types from petsc4py, you should do:

from petsc4py.PETSc cimport Vec,  PetscVec
from petsc4py.PETSc cimport Mat,  PetscMat
# etcetera

now you can "cdef"  either "PetscVec" (the low-level PETSc C type) and
"Vec" (the higher-level petsc4py Python type)

Or you can just

from petsc4py.PETSc cimport *

to get all definitions, or even

from petsc4py cimport PETSc

then use properly-namespaced definitions "PETSc.PetscVec" (C type) and
"PETSc.Vec" (Python type).


-- 
Lisandro Dalcin

Research Scientist
Computer, Electrical and Mathematical Sciences & Engineering (CEMSE)
Extreme Computing Research Center (ECRC)
King Abdullah University of Science and Technology (KAUST)
http://ecrc.kaust.edu.sa/

4700 King Abdullah University of Science and Technology
al-Khawarizmi Bldg (Bldg 1), Office # 0109
Thuwal 23955-6900, Kingdom of Saudi Arabia
http://www.kaust.edu.sa

Office Phone: +966 12 808-0459


Re: [petsc-users] Installing petsc4py

2017-10-18 Thread Lisandro Dalcin
Dear Nicholas,

The errors you reported where just outdated Fortran code and makefiles
in demo/, these are just examples. I pushed a fix to the maint branch
of the repository.

BTW, do you know that conda-forge provides prebuilt packages for petsc
and petsc4py? Just create a new conda environment, activate, and then
"conda install -c conda-forge petsc4py"

On 18 October 2017 at 07:23, Stegmeier, Nicholas
<nicholas.stegme...@sdstate.edu> wrote:
> Hello all,
>
>
> I was hoping to get some advice on installing petsc4py. Hopefully this is an
> appropriate question as it's my first time emailing the list.
>
>
> I failed numerous times installing petsc4py on my personal laptop, so I
> decided to use a virtual machine with Ubuntu to try installing from a clean
> slate.
>
>
> I followed the steps listed on this webpage:
> https://gist.github.com/mrosemeier/088115b2e34f319b913a
>
> Install openmpi
> Install anaconda-python
> Install petsc
> Install petsc4py
>
>
> After doing all this, I entered the petsc4py/demo folder and used make to
> see if all the tests passed. But several of the tests promptly failed. Below
> are some of the errors I received, but I also attached the full error
> output. Maybe I have forgotten to assign some environment variable?
>
> Thank you,
> Nicholas Stegmeier
>
> inside .bashrc:
> export LD_LIBRARY_PATH=/home/nick/local/openmpi-1.10.5/lib:$LD_LIBRARY_PATH
> export PATH=/home/nick/local/openmpi-1.10.5/bin:$PATH
> export PATH="/home/nick/local/anaconda2/bin:$PATH"
> export PETSC_DIR=/home/nick/local/petsc.git
> export PETSC_ARCH=arch-python-linux-i686
>
> nick@nick-VirtualBox:~/local/petsc.git/bin$ echo $PATH
> /home/nick/local/anaconda2/bin:/home/nick/local/openmpi-1.10.5/bin:/home/nick/local/anaconda2/bin:/home/nick/local/openmpi-1.10.5/bin:/home/nick/bin:/home/nick/.local/bin:/usr/local/sbin:/usr/local/bin:/usr/sbin:/usr/bin:/sbin:/bin:/usr/games:/usr/local/games:/snap/bin
>
> nick@nick-VirtualBox:~/local/petsc.git/bin$ echo $PATH
> /home/nick/local/anaconda2/bin:/home/nick/local/openmpi-1.10.5/bin:/home/nick/local/anaconda2/bin:/home/nick/local/openmpi-1.10.5/bin:/home/nick/bin:/home/nick/.local/bin:/usr/local/sbin:/usr/local/bin:/usr/sbin:/usr/bin:/sbin:/bin:/usr/games:/usr/local/games:/snap/bin
>
> sample error output:
>
> make -f makefile.petsc \
> PETSC_DIR=/home/nick/local/petsc.git
> PETSC_ARCH=arch-python-linux-i686
> make[2]: Entering directory '/home/nick/local/petsc4py.git/demo/perftest'
> cApp.f90
> make[2]: c: Command not found
> makefile.petsc:15: recipe for target 'driver.exe' failed
> make[2]: [driver.exe] Error 127 (ignored)
> /home/nick/local/openmpi-1.10.5/bin/mpicc -c -Wall -Wwrite-strings
> -Wno-strict-aliasing -Wno-unknown-pragmas -fstack-protector
> -fvisibility=hidden -g -O   -I/home/nick/local/petsc.git/include
> -I/home/nick/local/petsc.git/arch-python-linux-i686/include
> `pwd`/driver.c
> /home/nick/local/openmpi-1.10.5/bin/mpicc -Wall -Wwrite-strings
> -Wno-strict-aliasing -Wno-unknown-pragmas -fstack-protector
> -fvisibility=hidden -g -O  -o driver.exe
> -Wl,-rpath,/home/nick/local/petsc.git/arch-python-linux-i686/lib
> -L/home/nick/local/petsc.git/arch-python-linux-i686/lib -lpetsc -llapack
> -lblas -lm -ldl  driver.o App.o
> gcc: error: App.o: No such file or directory
> makefile.petsc:15: recipe for target 'driver.exe' failed
> make[2]: [driver.exe] Error 1 (ignored)
>
> ...
>
> ./Bratu2D.F90:20:0:
>
>  #include "petsc/finclude/petscdef.h"
>  ^
> Fatal Error: petsc/finclude/petscdef.h: No such file or directory
> compilation terminated.
> ./Bratu2D.F90:20:0:
>
>  #include "petsc/finclude/petscdef.h"
>  ^
> Fatal Error: petsc/finclude/petscdef.h: No such file or directory
> compilation terminated.
> error: Command "/usr/bin/gfortran -Wall -g -fno-second-underscore
> -DF2PY_REPORT_ON_ARRAY_COPY=1
> -I/home/nick/local/petsc.git/arch-python-linux-i686/include
> -I/home/nick/local/petsc.git/include
> -I/home/nick/local/anaconda2/lib/python2.7/site-packages/petsc4py/include
> -I. -Ibuild/src.linux-i686-2.7/.
> -I/home/nick/local/anaconda2/lib/python2.7/site-packages/numpy/core/include
> -I/home/nick/local/anaconda2/include/python2.7 -c -c ./Bratu2D.F90 -o
> build/temp.linux-i686-2.7/Bratu2D.o -Jbuild/temp.linux-i686-2.7/
> -Ibuild/temp.linux-i686-2.7/" failed with exit status 1
> makefile:20: recipe for target 'Bratu2D.so' failed
> make[1]: *** [Bratu2D.so] Error 1
> make[1]: Leaving directory '/home/nick/local/petsc4py.git/demo/wrap-f2py'
> makefile:3: recipe for target 'all' failed
> make: [all] Error 2 (ignored)
>
>
> Thank you,
> Nicholas Stegmeier



-- 
Lisandro Dalcin

Resea

Re: [petsc-users] Trouble installing petsc4py in Anaconda environment

2017-10-12 Thread Lisandro Dalcin
  PETSC_DIR:
> /home/login/anaconda2/envs/myenv/lib/python2.7/site-packages/petsc
> PETSC_ARCH:
> version:  3.8.0 release
> integer-size: 32-bit
> scalar-type:  real
> precision:double
> language: CONLY
> compiler: /vendor/sgi/mpt/2.14r19/bin/mpicc
> linker:   /vendor/sgi/mpt/2.14r19/bin/mpicc
> building 'PETSc' extension
> creating build/temp.linux-x86_64-2.7
> creating build/temp.linux-x86_64-2.7/src
> /vendor/sgi/mpt/2.14r19/bin/mpicc -pthread -B
> /home/login/anaconda2/envs/myenv/compiler_compat -fPIC -Wall -Wwrite-strings
> -Wno-strict-aliasing -Wno-unknown-pragmas -fstack-protector -g -O -fPIC
> -fno-strict-aliasing -g -O2 -DNDEBUG -g -fwrapv -O3 -Wall
> -DPETSC_DIR=/home/login/anaconda2/envs/myenv/lib/python2.7/site-packages/petsc
> -I/home/login/anaconda2/envs/myenv/lib/python2.7/site-packages/petsc/include
> -Isrc/include
> -I/home/login/anaconda2/envs/myenv/lib/python2.7/site-packages/numpy/core/include
> -I/home/login/anaconda2/envs/myenv/include/python2.7 -c src/PETSc.c -o
> build/temp.linux-x86_64-2.7/src/PETSc.o
> In file included from
> /home/login/anaconda2/envs/myenv/lib/python2.7/site-packages/numpy/core/include/numpy/ndarraytypes.h:1788:0,
>  from
> /home/login/anaconda2/envs/myenv/lib/python2.7/site-packages/numpy/core/include/numpy/ndarrayobject.h:18,
>  from
> /home/login/anaconda2/envs/myenv/lib/python2.7/site-packages/numpy/core/include/numpy/arrayobject.h:4,
>  from src/include/petsc4py/numpy.h:11,
>  from src/petsc4py.PETSc.c:519,
>  from src/PETSc.c:3:
>
> /home/login/anaconda2/envs/myenv/lib/python2.7/site-packages/numpy/core/include/numpy/npy_1_7_deprecated_api.h:15:2:
> warning: #warning "Using deprecated NumPy API, disable it by " "#defining
> NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION" [-Wcpp]
>  #warning "Using deprecated NumPy API, disable it by " \
>   ^
> /vendor/sgi/mpt/2.14r19/bin/mpicc -pthread -B
> /home/login/anaconda2/envs/myenv/compiler_compat -fPIC -Wall -Wwrite-strings
> -Wno-strict-aliasing -Wno-unknown-pragmas -fstack-protector -g -O -fPIC
> -fno-strict-aliasing -g -O2 -DNDEBUG -g -fwrapv -O3 -Wall
> -DPETSC_DIR=/home/login/anaconda2/envs/myenv/lib/python2.7/site-packages/petsc
> -I/home/login/anaconda2/envs/myenv/lib/python2.7/site-packages/petsc/include
> -Isrc/include
> -I/home/login/anaconda2/envs/myenv/lib/python2.7/site-packages/numpy/core/include
> -I/home/login/anaconda2/envs/myenv/include/python2.7 -c src/libpetsc4py.c -o
> build/temp.linux-x86_64-2.7/src/libpetsc4py.o
> /vendor/sgi/mpt/2.14r19/bin/mpicc -pthread -B
> /home/login/anaconda2/envs/myenv/compiler_compat -fPIC -Wall -Wwrite-strings
> -Wno-strict-aliasing -Wno-unknown-pragmas -fstack-protector -g -O gcc
> -pthread -shared -B /home/login/anaconda2/envs/myenv/compiler_compat
> -fno-strict-aliasing -g -O2 -DNDEBUG -g -fwrapv -O3 -Wall
> -L/home/login/anaconda2/envs/myenv/lib
> -Wl,-rpath=/home/login/anaconda2/envs/myenv/lib -Wl,--no-as-needed
> -Wl,--sysroot=/ build/temp.linux-x86_64-2.7/src/PETSc.o
> build/temp.linux-x86_64-2.7/src/libpetsc4py.o
> -L/home/login/anaconda2/envs/myenv/lib/python2.7/site-packages/petsc/lib
> -L/home/login/anaconda2/envs/myenv/lib
> -Wl,-rpath,/home/login/anaconda2/envs/myenv/lib/python2.7/site-packages/petsc/lib
> -Wl,-rpath,/home/login/anaconda2/envs/myenv/lib -lpetsc -lpython2.7 -o
> build/lib.linux-x86_64-2.7/petsc4py/lib/PETSc.so
> gcc: error: gcc: No such file or directory
> error: command '/vendor/sgi/mpt/2.14r19/bin/mpicc' failed with exit
> status 1
>
> 
> Command "/home/login/anaconda2/envs/myenv/bin/python -u -c "import
> setuptools,
> tokenize;__file__='/tmp/pip-build-v6lDIk/petsc4py/setup.py';f=getattr(tokenize,
> 'open', open)(__file__);code=f.read().replace('\r\n',
> '\n');f.close();exec(compile(code, __file__, 'exec'))" install --record
> /tmp/pip-d6NYW8-record/install-record.txt
> --single-version-externally-managed --compile" failed with error code 1 in
> /tmp/pip-build-v6lDIk/petsc4py/
>
>
> --
> =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
>
> Bill Jones   w.t.jo...@nasa.gov
> Mail Stop 128 Computational AeroSciences Branch
> 15 Langley Boulevard   Research Directorate
> NASA Langley Research Center   Building 1268, Room 1044
> Hampton, VA  23681-2199   Phone +1 757 864-5318
> Fax +1 757 864-8816
>  http://fun3d.larc.nasa.gov



-- 
Lisandro Dalcin

Research Scientist
Computer, Electrical and Mathematical Sciences & Engineering (CEMSE)
Extreme Computing Research Center (ECRC)
King Abdullah University of Science and Technology (KAUST)
http://ecrc.kaust.edu.sa/

4700 King Abdullah University of Science and Technology
al-Khawarizmi Bldg (Bldg 1), Office # 0109
Thuwal 23955-6900, Kingdom of Saudi Arabia
http://www.kaust.edu.sa

Office Phone: +966 12 808-0459


Re: [petsc-users] [PETSc] How to read a Vector and Matrix from an ASCII file into PETSc?

2017-09-11 Thread Lisandro Dalcin
On 11 September 2017 at 19:21, Ryan Morley <ryantmorley7...@gmail.com> wrote:
> Dear all,
>
> I am working on a program in Petsc that will read a vector and matrix from
> an ASCII file and store it into a vector and matrix object in Petsc so I can
> eventually solve. It needs to be an ASCII file because I will be outputting
> the vector and matrix to an ASCII file from code in Delphi.
>

Would it be OK to use an intermediate program to read the ASCII output
and dump in PETSc binary format? That would be rather trivial to do
with a Python script.


-- 
Lisandro Dalcin

Research Scientist
Computer, Electrical and Mathematical Sciences & Engineering (CEMSE)
Extreme Computing Research Center (ECRC)
King Abdullah University of Science and Technology (KAUST)
http://ecrc.kaust.edu.sa/

4700 King Abdullah University of Science and Technology
al-Khawarizmi Bldg (Bldg 1), Office # 0109
Thuwal 23955-6900, Kingdom of Saudi Arabia
http://www.kaust.edu.sa

Office Phone: +966 12 808-0459


Re: [petsc-users] TS for explicit equations

2017-07-07 Thread Lisandro Dalcin
On 6 July 2017 at 17:10, Alejandro D Otero <aot...@fi.uba.ar> wrote:
>
> # Set rhs function
> ts.setRHSFunction(self.evalRHSFunction, self.g)
> (where evalRHSFunction has the form f(self, ts, t, Vort, f))
>

Please double check that the RHS function "outputs" the result in the
passed-in vector "f", my guess is that you are just returning a new
vector that is eventually ignored.


-- 
Lisandro Dalcin

Research Scientist
Computer, Electrical and Mathematical Sciences & Engineering (CEMSE)
Extreme Computing Research Center (ECRC)
King Abdullah University of Science and Technology (KAUST)
http://ecrc.kaust.edu.sa/

4700 King Abdullah University of Science and Technology
al-Khawarizmi Bldg (Bldg 1), Office # 0109
Thuwal 23955-6900, Kingdom of Saudi Arabia
http://www.kaust.edu.sa

Office Phone: +966 12 808-0459


Re: [petsc-users] petsc4py and python's logging module

2017-05-29 Thread Lisandro Dalcin
On 29 May 2017 at 19:17, David Nolte <dno...@dim.uchile.cl> wrote:
> Dear all,
>
> is it possible to use python's logging module
> (https://docs.python.org/2/howto/logging.html) to handle PETSc output in
> python, such as the residuals during a KSP/SNES solve?
> I log my solver's activity to a file using the logging module, it would
> be great to include the PETSc output aswell.
>

Not sure if this is what you really want, but you could...

1) Use {ksp|snes}.setConvergenceHistory() before solve, then
{ksp|snes}.getConvergenceHistory() after solve, you will get arrays
with the residual history, then you can do whatever you want with
them.

2) Implement a KSP/SNES monitor in a Python function and call
{ksp|snes}.setMonitor(), then you can use python's logging inside your
monitor function.


-- 
Lisandro Dalcin

Research Scientist
Computer, Electrical and Mathematical Sciences & Engineering (CEMSE)
Extreme Computing Research Center (ECRC)
King Abdullah University of Science and Technology (KAUST)
http://ecrc.kaust.edu.sa/

4700 King Abdullah University of Science and Technology
al-Khawarizmi Bldg (Bldg 1), Office # 0109
Thuwal 23955-6900, Kingdom of Saudi Arabia
http://www.kaust.edu.sa

Office Phone: +966 12 808-0459


Re: [petsc-users] petsc4py bool type

2017-04-26 Thread Lisandro Dalcin
On 25 April 2017 at 19:48, Garth N. Wells <gn...@cam.ac.uk> wrote:
> I'm seeing some behaviour with bool types in petsc4py that I didn't
> expect. In the Python interface, returned Booleans have type ' 'int'>', where I expected them to have type ' '. Below
> program illustrates issue. Seems to be related to bint in cython. Am I
> doing something wrong?

Wow! Yes, you are right, this seems like a regression in Cython. In
the past, a  cast used to coerce values to the Python `bool`
type. Damn, these casts are everywhere, it might take a while to fix
all the code, unless I can find a hacky way to workaround the issue.


-- 
Lisandro Dalcin

Research Scientist
Computer, Electrical and Mathematical Sciences & Engineering (CEMSE)
Extreme Computing Research Center (ECRC)
King Abdullah University of Science and Technology (KAUST)
http://ecrc.kaust.edu.sa/

4700 King Abdullah University of Science and Technology
al-Khawarizmi Bldg (Bldg 1), Office # 0109
Thuwal 23955-6900, Kingdom of Saudi Arabia
http://www.kaust.edu.sa

Office Phone: +966 12 808-0459


Re: [petsc-users] petsc4py & GPU

2017-04-20 Thread Lisandro Dalcin
On 20 April 2017 at 17:09, Damian Kaliszan <dam...@man.poznan.pl> wrote:

> Thank you for reply:) Sorry for maybe stupid question in the scope of
> setting petsc(4py) options.
> Should the following calls (somewhere before creating matrix & vectors):
>
> PETSc.Options().setValue("ksp_view", "")
> PETSc.Options().setValue("log_view", "")
>

Unfortunately, no. There are a few options (-log_view ?) that you should
set before calling PetscInitialize() (which happens automatically at import
time), otherwise things do not work as expected. To pass things from the
command line and set them before PetscInitialize() the usual idiom is:

import sys, petsc4py
petsc4py.init(sys.argv)
from petsc4py import PETSc



-- 
Lisandro Dalcin

Research Scientist
Computer, Electrical and Mathematical Sciences & Engineering (CEMSE)
Extreme Computing Research Center (ECRC)
King Abdullah University of Science and Technology (KAUST)
http://ecrc.kaust.edu.sa/

4700 King Abdullah University of Science and Technology
al-Khawarizmi Bldg (Bldg 1), Office # 0109
Thuwal 23955-6900, Kingdom of Saudi Arabia
http://www.kaust.edu.sa

Office Phone: +966 12 808-0459


Re: [petsc-users] how to use petsc4py with mpi subcommunicators?

2017-04-13 Thread Lisandro Dalcin
On 11 April 2017 at 17:31, Rodrigo Felicio <rodrigo.feli...@iongeo.com>
wrote:

> Thanks, Jed, but using color == 0 lead to the same error msg. Is there no
> way to set PETSc.COMM_WORLD to a subcomm instead of MPI.COMM_WORLD in
> python?
>
>

You can do it, but it is a bit tricky, and related to the fact that
PETSC_COMM_WORLD has to be set to subcomm before calling PetscInitialize(),
which happens automatically the first time Python executes "from petsc4py
import PETSc". The way to do it would be the following:

import sys, petsc4py
petsc4py.init(sys.argv, comm=subcomm)
from petsc4py import PETSc



-- 
Lisandro Dalcin

Research Scientist
Computer, Electrical and Mathematical Sciences & Engineering (CEMSE)
Extreme Computing Research Center (ECRC)
King Abdullah University of Science and Technology (KAUST)
http://ecrc.kaust.edu.sa/

4700 King Abdullah University of Science and Technology
al-Khawarizmi Bldg (Bldg 1), Office # 0109
Thuwal 23955-6900, Kingdom of Saudi Arabia
http://www.kaust.edu.sa

Office Phone: +966 12 808-0459


Re: [petsc-users] How to use f2py on a PETSc/SLEPc-based fortran code

2017-03-25 Thread Lisandro Dalcin
On 22 March 2017 at 20:29, Barry Smith <bsm...@mcs.anl.gov> wrote:

>
>   Lisandro,
>
> We've had a couple questions similar to this with f2py; is there a way
> we could add to the PETSc/SLEPc makefile rules something to allow people to
> trivially use f2py without having to make their own (often incorrect)
> manual command lines?
>
>Thanks
>
>
Barry, it is quite hard and hacky to get f2py working in the general case.
I think the email from Gaetan in this thread proves my point.

IMHO, it is easier to write a small Fortran source exposing the API to call
using ISO_C_BINDINGS, then wrap that code with the more traditional C-based
"static" tools (SWIG, Cython) or even "dynamically" with ctypes or cffi
(which use dlopen'ing).



-- 
Lisandro Dalcin

Research Scientist
Computer, Electrical and Mathematical Sciences & Engineering (CEMSE)
Extreme Computing Research Center (ECRC)
King Abdullah University of Science and Technology (KAUST)
http://ecrc.kaust.edu.sa/

4700 King Abdullah University of Science and Technology
al-Khawarizmi Bldg (Bldg 1), Office # 0109
Thuwal 23955-6900, Kingdom of Saudi Arabia
http://www.kaust.edu.sa

Office Phone: +966 12 808-0459


Re: [petsc-users] FFT using Petsc4py

2016-09-30 Thread Lisandro Dalcin
On 27 September 2016 at 21:05, Amit Itagi <gcf...@gmail.com> wrote:

> Hello,
>
> I am looking at the Petsc FFT interfaces. I was wondering if a parallel
> FFT can be performed within a Petsc4Py code. If not, what would be the best
> way to use the Petsc interfaces for FFT from Petsc4Py ?
>
>
It should work out of the box by using mat.setType(Mat.Type.FFTW) before
setup of your matrix.


-- 
Lisandro Dalcin

Research Scientist
Computer, Electrical and Mathematical Sciences & Engineering (CEMSE)
Extreme Computing Research Center (ECRC)
King Abdullah University of Science and Technology (KAUST)
http://ecrc.kaust.edu.sa/

4700 King Abdullah University of Science and Technology
al-Khawarizmi Bldg (Bldg 1), Office # 0109
Thuwal 23955-6900, Kingdom of Saudi Arabia
http://www.kaust.edu.sa

Office Phone: +966 12 808-0459


Re: [petsc-users] petsc4py with complex numbers?

2016-09-26 Thread Lisandro Dalcin
Are you sure you built petsc4py with a PETSc build with complex scalars?
Whats the output of  "print(PETSc.ScalarType)" ?


On 26 September 2016 at 07:47, Aurelien Ponte <aurelien.po...@ifremer.fr>
wrote:

> Hi,
>
> I am trying to solve a linear problem whose operator has complex
> coefficients via petsc4py
> but keep running into the following error:
>
>   File 
> "/Users/aponte/Current_projects/people/kraig_marine/wd_response/solver/set_L.py",
> line 52, in set_L
> L.setValueStencil(row, col, value)
>   File "PETSc/Mat.pyx", line 882, in petsc4py.PETSc.Mat.setValueStencil
> (src/petsc4py.PETSc.c:124785)
>   File "PETSc/petscmat.pxi", line 1017, in petsc4py.PETSc.matsetvaluestencil
> (src/petsc4py.PETSc.c:31469)
>   File "PETSc/arraynpy.pxi", line 140, in petsc4py.PETSc.iarray_s
> (src/petsc4py.PETSc.c:8811)
>   File "PETSc/arraynpy.pxi", line 121, in petsc4py.PETSc.iarray
> (src/petsc4py.PETSc.c:8542)
> TypeError: can't convert complex to float
>
> The code looks like:
>
> for j in range(ys, ye):
> for i in range(xs, xe):
> row.index = (i, j, kx)
> row.field = 0
> col.index = (i, j, kx) col.field=0
> L.setValueStencil(row, col, 1j)
>
> Any idea about what am I doing wrong?
>
> cheers
>
> aurelien
>
>
>
> --
> Aurélien Ponte
> Tel: (+33) 2 98 22 40 73
> Fax: (+33) 2 98 22 44 96
> UMR 6523, IFREMER
> ZI de la Pointe du Diable
> CS 10070
> 29280 Plouzané
>
>


-- 
Lisandro Dalcin

Research Scientist
Computer, Electrical and Mathematical Sciences & Engineering (CEMSE)
Extreme Computing Research Center (ECRC)
King Abdullah University of Science and Technology (KAUST)
http://ecrc.kaust.edu.sa/

4700 King Abdullah University of Science and Technology
al-Khawarizmi Bldg (Bldg 1), Office # 0109
Thuwal 23955-6900, Kingdom of Saudi Arabia
http://www.kaust.edu.sa

Office Phone: +966 12 808-0459


Re: [petsc-users] [petsc4py] a problem with computeRHSFunctionLinear interface?

2016-08-21 Thread Lisandro Dalcin
On 21 August 2016 at 13:01, Francesco Caimmi <francesco.cai...@polimi.it>
wrote:

> I will also take the chance to ask a related question: in the original C
> code,
> lines 213 -214 read
> ierr = TSGetSNES(ts,);CHKERRQ(ierr);
> ierr=SNESSetJacobian(snes,NULL,NULL,SNESComputeJacobianDefault,
> NULL);CHKERRQ(ierr);
>
> How do I translate the last line in Python? I was not able to find the
> equivalent of SNESComputeJacobianDefault.
>

That one is not directly available, however you can pass -snes_fd in the
command line, or alternatively

opts = PETSc.Options()
opts['snes_fd'] = 1
snes.setFromOptions()


-- 
Lisandro Dalcin

Research Scientist
Computer, Electrical and Mathematical Sciences & Engineering (CEMSE)
Extreme Computing Research Center (ECRC)
King Abdullah University of Science and Technology (KAUST)
http://ecrc.kaust.edu.sa/

4700 King Abdullah University of Science and Technology
al-Khawarizmi Bldg (Bldg 1), Office # 0109
Thuwal 23955-6900, Kingdom of Saudi Arabia
http://www.kaust.edu.sa

Office Phone: +966 12 808-0459


Re: [petsc-users] [petsc4py] a problem with computeRHSFunctionLinear interface?

2016-08-18 Thread Lisandro Dalcin
Dear Francesco, sorry for the late answer, I missed your email.

You have to do it this way:

ts.setRHSFunction(PETSc.TS.computeRHSFunctionLinear)
ts.setRHSJacobian(PETSc.TS.computeRHSJacobianConstant, J=A, P=A)

I.e, you have to set unbound methods, not instance methods as you did
in your original code. Additionally, do not pass "args" nor "kargs".


On 11 August 2016 at 10:36, Francesco Caimmi <francesco.cai...@polimi.it> wrote:
> Dear all,
>
> I was trying to reproduce /ts/examples/tutorials/ex4.c in python to learn how
> to use TS solvers; the example uses the function TSComputeRHSFunctionLinear.
> However I get an error when running my code (attached in case you want to look
> at it), when I call ts.solve.
>
> Here is the trace:
> [fcaimmi@Wotan 2645] > ./ts_ex4.py
> Solving a linear TS problem, number of processors =  1
> Timestep  0 : time = 0.0  2-norm error = 1.14956855594e-08  max norm error = 0
> Traceback (most recent call last):
>   File "./ts_ex4.py", line 473, in 
> main(m = m, debug = debug)
>   File "./ts_ex4.py", line 340, in main
> ts.solve(u)
>   File "PETSc/TS.pyx", line 568, in petsc4py.PETSc.TS.solve
> (src/petsc4py.PETSc.c:188927)
>   File "PETSc/petscts.pxi", line 221, in petsc4py.PETSc.TS_RHSFunction
> (src/petsc4py.PETSc.c:35490)
>   File "PETSc/TS.pyx", line 189, in petsc4py.PETSc.TS.computeRHSFunctionLinear
> (src/petsc4py.PETSc.c:181611)
> TypeError: computeRHSFunctionLinear() takes exactly 3 positional arguments (5
> given)
>
> I cannot understand if there is a problem with my code or if the problem is in
> computeRHSFunctionLinear interface.
>
>  I checked https://bitbucket.org/petsc/petsc4py/ and the interface to
> computeRHSFunctionLinear  has three arguments, however I am not that much into
> petsc4py to tell how it gets called.
>
> I am on Petsc Release Version 3.7.3
>
> Thank you for your time.
>
> Best,
> --
> Francesco Caimmi
>
> Laboratorio di Ingegneria dei Polimeri
> http://www.chem.polimi.it/polyenglab/
>
> Politecnico di Milano - Dipartimento di Chimica,
> Materiali e Ingegneria Chimica “Giulio Natta”
>
> P.zza Leonardo da Vinci, 32
> I-20133 Milano
> Tel. +39.02.2399.4711
> Fax +39.02.7063.8173
>
> francesco.cai...@polimi.it
> Skype: fmglcaimmi (please arrange meetings by e-mail)



-- 
Lisandro Dalcin

Research Scientist
Computer, Electrical and Mathematical Sciences & Engineering (CEMSE)
Extreme Computing Research Center (ECRC)
King Abdullah University of Science and Technology (KAUST)
http://ecrc.kaust.edu.sa/

4700 King Abdullah University of Science and Technology
al-Khawarizmi Bldg (Bldg 1), Office # 0109
Thuwal 23955-6900, Kingdom of Saudi Arabia
http://www.kaust.edu.sa

Office Phone: +966 12 808-0459


Re: [petsc-users] petsc4py - Change default line search for SNES Newton Line Search

2016-07-12 Thread Lisandro Dalcin
On 12 July 2016 at 16:37, Hensgens, Marius
<marius.hensg...@rwth-aachen.de> wrote:
> how can I change the used line search method after setting the SNES Type to
> 'newtonls' using petsc4py ?
>

Right now, you can either use the command line or programatically
insert an option in the database

opts = PETSc.Options()
opts['snes_linesearch_type'] = lstype
...
snes.setFromOptions()

> In the official PETSc documentation there is a function called
> SNESLineSearchSetType, however in petsc4py I can't find an equivalent
> function.

The SNESLineSearch type and related routines are not wrapped yet.


-- 
Lisandro Dalcin

Research Scientist
Computer, Electrical and Mathematical Sciences & Engineering (CEMSE)
Extreme Computing Research Center (ECRC)
King Abdullah University of Science and Technology (KAUST)
http://ecrc.kaust.edu.sa/

4700 King Abdullah University of Science and Technology
al-Khawarizmi Bldg (Bldg 1), Office # 0109
Thuwal 23955-6900, Kingdom of Saudi Arabia
http://www.kaust.edu.sa

Office Phone: +966 12 808-0459


Re: [petsc-users] How to determine the type of SNESLineSearch?

2016-07-10 Thread Lisandro Dalcin
PetscBool match;
ierr = PetscObjectTypeCompare((PetscObject)linesearch,
SNESLINESEARCHBASIC,);CHKERRQ(ierr);
if (!match) {
...
}

On 6 July 2016 at 00:07, Gautam Bisht <gbi...@lbl.gov> wrote:
> Hi PETSc,
>
> After SNESSolve converges, I want to perform few additional operations only
> when SNESLineSearchType is not SNESLINESEARCHBASIC. But, there is no
> SNESLineSearchGetType routine. Any idea on how I can determine the type of
> LineSearch set by a user using command line option?
>
> Thanks,
> -Gautam.



-- 
Lisandro Dalcin

Research Scientist
Computer, Electrical and Mathematical Sciences & Engineering (CEMSE)
Extreme Computing Research Center (ECRC)
King Abdullah University of Science and Technology (KAUST)
http://ecrc.kaust.edu.sa/

4700 King Abdullah University of Science and Technology
al-Khawarizmi Bldg (Bldg 1), Office # 0109
Thuwal 23955-6900, Kingdom of Saudi Arabia
http://www.kaust.edu.sa

Office Phone: +966 12 808-0459


Re: [petsc-users] Petsc4py DMREDUNDANT or equivalent functionality

2016-04-15 Thread Lisandro Dalcin
On 15 April 2016 at 20:17, Ben Fogelson <benfogel...@gmail.com> wrote:
> As best I can tell, Petsc4py doesn't have a method like
> DM().createRedundant()--perhaps it this DM just isn't implemented?

Yes, that is the case.

> Can
> anyone recommend an alternate approach to getting the same functionality?

I guess the best way to go would be to implement the missing DM.


-- 
Lisandro Dalcin

Research Scientist
Computer, Electrical and Mathematical Sciences & Engineering (CEMSE)
Extreme Computing Research Center (ECRC)
King Abdullah University of Science and Technology (KAUST)
http://ecrc.kaust.edu.sa/

4700 King Abdullah University of Science and Technology
al-Khawarizmi Bldg (Bldg 1), Office # 4332
Thuwal 23955-6900, Kingdom of Saudi Arabia
http://www.kaust.edu.sa

Office Phone: +966 12 808-0459


Re: [petsc-users] PETSc release soon, request for input on needed fixes or enhancements

2016-02-29 Thread Lisandro Dalcin
On 27 February 2016 at 23:36, Barry Smith <bsm...@mcs.anl.gov> wrote:
>   We are planning the PETSc release 3.7 shortly. If you know of any bugs that 
> need to be fixed or enhancements added before the release please let us know.

Barry, long time ago I reported a regression in PetIGA related to the
new Vec assembly routines. I guess Jed didn't have a chance to look at
it. Until that happens, I propose to not use the new routines as the
default ones.


-- 
Lisandro Dalcin

Research Scientist
Computer, Electrical and Mathematical Sciences & Engineering (CEMSE)
Extreme Computing Research Center (ECRC)
King Abdullah University of Science and Technology (KAUST)
http://ecrc.kaust.edu.sa/

4700 King Abdullah University of Science and Technology
al-Khawarizmi Bldg (Bldg 1), Office # 4332
Thuwal 23955-6900, Kingdom of Saudi Arabia
http://www.kaust.edu.sa

Office Phone: +966 12 808-0459


Re: [petsc-users] petsc4py installation error with PETSc v 3.6.1

2015-07-25 Thread Lisandro Dalcin
On 25 July 2015 at 00:12, Justin Chang jychan...@gmail.com wrote:
 ImportError:
 /usr/local/lib/python2.7/dist-packages/petsc4py/lib/arch-linux2-c-opt/PETSc.so:
 undefined symbol: DMShellSetCreateSubDM
 make: *** [unit_sequential] Error 1


You need to git pull in branch master of petsc-dev and rebuild.


-- 
Lisandro Dalcin

Research Scientist
Computer, Electrical and Mathematical Sciences  Engineering (CEMSE)
Numerical Porous Media Center (NumPor)
King Abdullah University of Science and Technology (KAUST)
http://numpor.kaust.edu.sa/

4700 King Abdullah University of Science and Technology
al-Khawarizmi Bldg (Bldg 1), Office # 4332
Thuwal 23955-6900, Kingdom of Saudi Arabia
http://www.kaust.edu.sa

Office Phone: +966 12 808-0459


Re: [petsc-users] petsc4py Build Problem

2015-06-24 Thread Lisandro Dalcin
On 23 June 2015 at 12:15, Mikhail Khodak mkho...@princeton.edu wrote:
 Thank you for reviewing this. The errors for the command 'pip2 install .'
 inside the petsc4py directory are attached, along with the configure/make
 logs for the PETSc build.
 Best,
 Mikhail Khodak


Not sure what's going on.

Do you have PETSC_DIR and PETSC_ARCH in your environment?

Does it work if you just do python setup.py build ?


-- 
Lisandro Dalcin

Research Scientist
Computer, Electrical and Mathematical Sciences  Engineering (CEMSE)
Numerical Porous Media Center (NumPor)
King Abdullah University of Science and Technology (KAUST)
http://numpor.kaust.edu.sa/

4700 King Abdullah University of Science and Technology
al-Khawarizmi Bldg (Bldg 1), Office # 4332
Thuwal 23955-6900, Kingdom of Saudi Arabia
http://www.kaust.edu.sa

Office Phone: +966 12 808-0459


Re: [petsc-users] petsc4py Build Problem

2015-06-22 Thread Lisandro Dalcin
Dear Mikhail, Satish helped me and we managed to get petsc4py working
with static builds and Cygwin's Open MPI. Please git pull and try
again.

On 11 June 2015 at 14:13, Mikhail Khodak mkho...@princeton.edu wrote:
 I have attached the log files.
 I tried using Open MPI before but received the same result. I will try a
 static build later and otherwise will get a VM.
 Thanks for looking at this,
 Mikhail


 On Thu, Jun 11, 2015 at 10:48 AM, Satish Balay ba...@mcs.anl.gov wrote:

 Also - should mention:

 - cygwin has OpenMPI pacakged - so you could also try that.
 - most of my builds are static [so never tried dll builds]

 Satish

 On Thu, 11 Jun 2015, Satish Balay wrote:

  If you get errors when running basic petsc examples - send us the
  relavant petsc logs [configure.log, make.log, test.log etc..]
 
  Also - note that mpich is not supported on cygwin/windows [it
  generally works for us - when we try the --download-mpich build].
 
  Unless you really need a cywin build/use of petsc4py - it might be
  easier to install a linux VM - and use PETSc/petsc4py on it.
 
  Satish
 
  On Thu, 11 Jun 2015, Mikhail Khodak wrote:
 
   I added these lines to the petsc4py test files (specifically
   test_comm.py)
   but the error remains the same. However, I have done the standard
   'which
   mpiexec', 'which mpicc', 'which mpirun' and they are all in the same
   folder. In fact it is the only MPI installed.
   The reason I thought it might be a PETSc build problem is because one
   of
   the PETSc 'make test' tests (ex5f) fails with the same error, even
   though
   the 'make streams' test works fine with MPI processes.
   Thanks,
   Mikhail
  
   On Thu, Jun 11, 2015 at 4:21 AM, Matthew Knepley knep...@gmail.com
   wrote:
  
On Thu, Jun 11, 2015 at 5:07 AM, Mikhail Khodak
mkho...@princeton.edu
wrote:
   
Thank you for your help - the install seems to work, apart from
routines
requiring MPI, which fail due to the Attempting to use an MPI
routine
before initializing MPI error. This seems to be an error in the
PETSc
build itself.
   
   
The MPI routines will not work until after
   
import petsc4py, sys
petsc4py.init(sys.argv)
from petsc4py import PETSc
   
If they fail after this, it is usually a mismatch between the
mpiexec
being used and the MPI
libraries being linked.
   
  Thanks,
   
Matt
   
Thanks again,
Mikhail Khodak
   
On Mon, Jun 8, 2015 at 5:11 AM, Lisandro Dalcin dalc...@gmail.com
wrote:
   
On 8 June 2015 at 02:50, Mikhail Khodak mkho...@princeton.edu
wrote:
 Hello, I am trying to build petsc4py-3.5.1 using Cygwin on
 64-bit
Windows 7.
 I have built PETSc 3.5.4 with shared and dynamic libraries using
 mpich2-1.2.1 and successfully ran the installation tests. I am
 using
Python
 2.7 and NumPy 1.9.2 and have installed mpi4py. However, when I
 attempt
to
 install petsc4py (both with pip and distutils) I get a mpicc
 compiler
error
 due to undefined references/symbols. I have attached the output
 of
running

 pip install petsc petsc4py --allow-external petsc
 --allow-external
petsc4py

   
I've never ever built or test petsc4py under Cygwin. The errors
you
see are expected.
Perhaps you can manually workaround the issues following the
following
steps:
   
1) Download the petsc4py tarball and unpack it.
2) Open the file src/libpetsc4py/libpetsc4py.h, add remove all
occurences of DL_IMPORT, i.e, replace DL_IMPORT(XYZ) for just XYZ
3) Use pip again:
   
pip install petsc
pip install .
   
The last line assumes your current working directory is the one
having
petsc4py's setup.py
   
Finally, I do not guarantee this will work. I'm just guessing,
petsc4py never explicitly supported Windows and/or Cygwin.
   
   
--
Lisandro Dalcin

Research Scientist
Computer, Electrical and Mathematical Sciences  Engineering
(CEMSE)
Numerical Porous Media Center (NumPor)
King Abdullah University of Science and Technology (KAUST)
http://numpor.kaust.edu.sa/
   
4700 King Abdullah University of Science and Technology
al-Khawarizmi Bldg (Bldg 1), Office # 4332
Thuwal 23955-6900, Kingdom of Saudi Arabia
http://www.kaust.edu.sa
   
Office Phone: +966 12 808-0459
   
   
   
   
   
--
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
   
  
 
 



-- 
Lisandro Dalcin

Research Scientist
Computer, Electrical and Mathematical Sciences  Engineering (CEMSE)
Numerical Porous Media Center (NumPor)
King Abdullah University of Science and Technology (KAUST)
http://numpor.kaust.edu.sa/

4700 King Abdullah University of Science and Technology
al-Khawarizmi

Re: [petsc-users] petsc4py Build Problem

2015-06-08 Thread Lisandro Dalcin
On 8 June 2015 at 02:50, Mikhail Khodak mkho...@princeton.edu wrote:
 Hello, I am trying to build petsc4py-3.5.1 using Cygwin on 64-bit Windows 7.
 I have built PETSc 3.5.4 with shared and dynamic libraries using
 mpich2-1.2.1 and successfully ran the installation tests. I am using Python
 2.7 and NumPy 1.9.2 and have installed mpi4py. However, when I attempt to
 install petsc4py (both with pip and distutils) I get a mpicc compiler error
 due to undefined references/symbols. I have attached the output of running

 pip install petsc petsc4py --allow-external petsc --allow-external petsc4py


I've never ever built or test petsc4py under Cygwin. The errors you
see are expected.
Perhaps you can manually workaround the issues following the following steps:

1) Download the petsc4py tarball and unpack it.
2) Open the file src/libpetsc4py/libpetsc4py.h, add remove all
occurences of DL_IMPORT, i.e, replace DL_IMPORT(XYZ) for just XYZ
3) Use pip again:

pip install petsc
pip install .

The last line assumes your current working directory is the one having
petsc4py's setup.py

Finally, I do not guarantee this will work. I'm just guessing,
petsc4py never explicitly supported Windows and/or Cygwin.


-- 
Lisandro Dalcin

Research Scientist
Computer, Electrical and Mathematical Sciences  Engineering (CEMSE)
Numerical Porous Media Center (NumPor)
King Abdullah University of Science and Technology (KAUST)
http://numpor.kaust.edu.sa/

4700 King Abdullah University of Science and Technology
al-Khawarizmi Bldg (Bldg 1), Office # 4332
Thuwal 23955-6900, Kingdom of Saudi Arabia
http://www.kaust.edu.sa

Office Phone: +966 12 808-0459


Re: [petsc-users] Create identity IS

2015-04-20 Thread Lisandro Dalcin
On 20 April 2015 at 14:35, Florian Lindner mailingli...@xgm.de wrote:
 Hello,

 I use an index set for a row mapping but don't want to use one for the column 
 mapping.

 I try to create an idenitity index set to supply to 
 ISLocalToGlobalMappingCreateIS and MatSetLocalToGlobalMapping:

   IS is;
   ISCreate(PETSC_COMM_WORLD, is);

At this point the index set does not have the type set.

   ISSetIdentity(is);

 but I get an segmentation fault at the last line.


Bad error checking, you should get an error saying that the type is not yet set.

 What is the best way to create an identity index set?


I would use ISCreateStride(), and perhaps next ISToGeneral() if you
really want a ISGENERAL type.


-- 
Lisandro Dalcin

Research Scientist
Computer, Electrical and Mathematical Sciences  Engineering (CEMSE)
Numerical Porous Media Center (NumPor)
King Abdullah University of Science and Technology (KAUST)
http://numpor.kaust.edu.sa/

4700 King Abdullah University of Science and Technology
al-Khawarizmi Bldg (Bldg 1), Office # 4332
Thuwal 23955-6900, Kingdom of Saudi Arabia
http://www.kaust.edu.sa

Office Phone: +966 12 808-0459


Re: [petsc-users] TimeStepper norm problems. EMIL Please read this

2015-04-16 Thread Lisandro Dalcin
On 16 April 2015 at 04:51, Emil Constantinescu emcon...@mcs.anl.gov wrote:

 On 3/24/15 5:31 AM, Lisandro Dalcin wrote:

 Emil, is there any chance you can try my code with your problem? I
 really need some feedback to push this to PETSc, otherwise


 Hi Lisandro - we checked ts_alpha_adapt and we tested it on a small system
 (mildly stiff van der Pol ODE). I enclosed a Figure generated by Debo that
 compares the error at the final time against ATOL - there alpha is the
 original one (with original adaptor).

Original adaptor? Do you mean TSAlphaSetAdapt(ts,TSAlphaAdaptDefault,NULL)?

In such a case, no surprises, TSAlphaAdaptDefault() is quite naive, it
does not estimate the LTE, it is actually based in the change of the
solution.

Can you send/point me your code to take a look?

 For this problem your adaptor works
 better [ error should be close to ATOL].

 So from my perspective, your
 adaptor works well. However, we need to keep the old too so that we have a
 one step adaptor available ... since alpha is a one step method.


Sorry, I'm not following you. Generalized alpha is not truly a
one-step method, because it retains Xdot from the previous time step.
In the initial step, other codes out there assume that Xdot_prev=0.

My alpha1 implementation has a trick for the initial time step
(using the same ideas, i.e. a BD formula) to both estimate Xdot_prev
and the error. However, this is at the cost of performing two
nonlinear solves in the initial step (IOW, two internal timesteps with
backward Euler and dt/2). Please note this trick should be valid for
any second-order method, at least to estimate the initial LTE in the
first step.

BTW, Have you compared my -ts_type tsalpha1 -ts_alpha_adapt with
PETSc's -ts_type cn -ts_theta adapt?

-- 
Lisandro Dalcin

Research Scientist
Computer, Electrical and Mathematical Sciences  Engineering (CEMSE)
Numerical Porous Media Center (NumPor)
King Abdullah University of Science and Technology (KAUST)
http://numpor.kaust.edu.sa/

4700 King Abdullah University of Science and Technology
al-Khawarizmi Bldg (Bldg 1), Office # 4332
Thuwal 23955-6900, Kingdom of Saudi Arabia
http://www.kaust.edu.sa

Office Phone: +966 12 808-0459


Re: [petsc-users] TimeStepper norm problems. EMIL Please read this

2015-04-16 Thread Lisandro Dalcin
On 16 April 2015 at 16:44, Emil Constantinescu emcon...@mcs.anl.gov wrote:
 On 4/16/15 2:13 AM, Lisandro Dalcin wrote:

 On 16 April 2015 at 04:51, Emil Constantinescu emcon...@mcs.anl.gov
 wrote:


 On 3/24/15 5:31 AM, Lisandro Dalcin wrote:


 Emil, is there any chance you can try my code with your problem? I
 really need some feedback to push this to PETSc, otherwise



 Hi Lisandro - we checked ts_alpha_adapt and we tested it on a small
 system
 (mildly stiff van der Pol ODE). I enclosed a Figure generated by Debo
 that
 compares the error at the final time against ATOL - there alpha is the
 original one (with original adaptor).


 Original adaptor? Do you mean
 TSAlphaSetAdapt(ts,TSAlphaAdaptDefault,NULL)?

 In such a case, no surprises, TSAlphaAdaptDefault() is quite naive, it
 does not estimate the LTE, it is actually based in the change of the
 solution.


 Yups, that's the one;

I was the author of TSAlphaAdaptDefault (the whole TSALPHA, in fact).
I coded that beast out of desperation long time ago. This thing is
not even implemented withing the current TSAdapt framework!. I think
at least we should re-implement this beast within TSAdapt (I mean,
similar as in -ts_theta_adapt).

 but that still needs to be available unless a better
 one-step one is implemented.

I think you have that feeling simply because TSAlphaAdaptDefault()
seems to work, but if you look carefully at the code, it does not make
sense from the point of view of a LTE-based theory.

 OK, Then let's define a poor's man one-step adapt that at least have
some sense. How would you implement adaptivity for backward Euler or
the midpoint rule? Something based exclusively in some l2/inf norm of
Xdot? IOW, we define LTE = X^{n+1} - X^n of O(\delta_t) and use the
usual formula with exponent order=1 ?

 Note however that even estimating exactly LTE
 is not a guarantee that the error will be within ATOL.


Well, ATOL is to keep the LTE under control, of course there are no
guarantees about the global error at the end step. That's what you
meant?


-- 
Lisandro Dalcin

Research Scientist
Computer, Electrical and Mathematical Sciences  Engineering (CEMSE)
Numerical Porous Media Center (NumPor)
King Abdullah University of Science and Technology (KAUST)
http://numpor.kaust.edu.sa/

4700 King Abdullah University of Science and Technology
al-Khawarizmi Bldg (Bldg 1), Office # 4332
Thuwal 23955-6900, Kingdom of Saudi Arabia
http://www.kaust.edu.sa

Office Phone: +966 12 808-0459


Re: [petsc-users] petsc4py and mumps (or other direct sparse solver)

2015-04-15 Thread Lisandro Dalcin
On 15 April 2015 at 18:15, Hong hzh...@mcs.anl.gov wrote:
 Benoit :

 mumps options can be called from PETSc C-code as
 ierr = PetscOptionsInsertString(-mat_mumps_icntl_13 1 -mat_mumps_icntl_24 1
 -mat_mumps_cntl_3 1e-12);CHKERRQ(ierr);

 or see petsc/src/ksp/ksp/examples/tutorials/ex52.c

 I do not know how to translate these calls into python API though.


opts = PETSc.Options() # dict-like
opts[mat_mumps_icntl_13] = 1
opts[mat_mumps_icntl_24] = 1
opts[mat_mumps_cntl_3] = 1e-12

Or also:

opts = PETSc.Options()
opts.prefixPush(mat_mumps_)
opts[icntl_13] = 1
opts[icntl_24] = 1
opts[cntl_3] = 1e-12
opts.prefixPop()

Or also:

opts = PETSc.Options(mat_mumps_)
opts[icntl_13] = 1
opts[icntl_24] = 1
opts[cntl_3] = 1e-12


-- 
Lisandro Dalcin

Research Scientist
Computer, Electrical and Mathematical Sciences  Engineering (CEMSE)
Numerical Porous Media Center (NumPor)
King Abdullah University of Science and Technology (KAUST)
http://numpor.kaust.edu.sa/

4700 King Abdullah University of Science and Technology
al-Khawarizmi Bldg (Bldg 1), Office # 4332
Thuwal 23955-6900, Kingdom of Saudi Arabia
http://www.kaust.edu.sa

Office Phone: +966 12 808-0459


Re: [petsc-users] running Petsc programs without a connection.

2015-04-01 Thread Lisandro Dalcin
On 26 March 2015 at 10:15, Sanjay Kharche
sanjay.khar...@manchester.ac.uk wrote:
 Is there a way of configuring or otherwise to get Petsc programs to run when 
 I do not have an internet connection?

Make sure you have localhost defined in /etc/hosts, I'm had the same
issue in Fedora with fresh installs.

$ cat /etc/hosts
127.0.0.1   localhost localhost.localdomain localhost4
localhost4.localdomain4 kw2060 kw2060.kaust.edu.sa
::1 localhost localhost.localdomain localhost6 localhost6.localdomain6



-- 
Lisandro Dalcin

Research Scientist
Computer, Electrical and Mathematical Sciences  Engineering (CEMSE)
Numerical Porous Media Center (NumPor)
King Abdullah University of Science and Technology (KAUST)
http://numpor.kaust.edu.sa/

4700 King Abdullah University of Science and Technology
al-Khawarizmi Bldg (Bldg 1), Office # 4332
Thuwal 23955-6900, Kingdom of Saudi Arabia
http://www.kaust.edu.sa

Office Phone: +966 12 808-0459


Re: [petsc-users] petsc4py: Passing command line switches

2015-03-26 Thread Lisandro Dalcin
On 25 March 2015 at 18:34, Matthew Knepley knep...@gmail.com wrote:
 That file is obviously old and was removed.


After installing petsc4py, you can do make docs in the top level
source tree, and the docs/ directory will be populated. For this to
work, you need sphinx and epydoc installed.

 Lisandro, are we generating documentation somewhere?

Two places:

+ http://petsc4py.readthedocs.org/ - This does not contain the
epydoc-generated API reference. I'm not sure if setting up readthedocs
to use epydoc is ever possible.

+ https://pythonhosted.org/petsc4py/,  - It is for the last release,
not the in-development version.



-- 
Lisandro Dalcin

Research Scientist
Computer, Electrical and Mathematical Sciences  Engineering (CEMSE)
Numerical Porous Media Center (NumPor)
King Abdullah University of Science and Technology (KAUST)
http://numpor.kaust.edu.sa/

4700 King Abdullah University of Science and Technology
al-Khawarizmi Bldg (Bldg 1), Office # 4332
Thuwal 23955-6900, Kingdom of Saudi Arabia
http://www.kaust.edu.sa

Office Phone: +966 12 808-0459


Re: [petsc-users] TimeStepper norm problems. EMIL Please read this

2015-03-24 Thread Lisandro Dalcin
On 24 March 2015 at 01:34, Jed Brown j...@jedbrown.org wrote:
 Lisandro Dalcin dalc...@gmail.com writes:

 On 21 March 2015 at 17:32, Emil Constantinescu emcon...@mcs.anl.gov wrote:
 When -ts_theta_adapt is used, then it detects the instability as an error
 and reduces the step by a lot! wlte=1.24e+03 which means that the reduction
 should be severe but the controller tries 0.1*dt and that seems to pass but
 it jig-saws (take a look at the next attempted step), which means that it
 is likely unstable.

 I think -ts_theta_adapt is seriously broken, I cannot make sense of
 the way the error estimator is computed.

 This was Shri's implementation.  TSEvaluateStep_Theta looks wrong to me,
 both using U as an input and assuming that Xdot is set.  I don't even
 know what the intended math is for a first-order embedded estimate in
 midpoint (theta=0.5).

Pease note that -ts_adapt_tetha errors for midpoint. You actually need
-ts_theta_endpoint to use adaptivity. But it is there, in the endpoint
version, that I cannot make sense of the implementation. The
trapezoidal rule is second-order, thus the LTE is O(dt^3). The
computed error estimator should be O(dt^2), but the current code seems
to compute an estimator O(dt), and that's the root of all evil.


-- 
Lisandro Dalcin

Research Scientist
Computer, Electrical and Mathematical Sciences  Engineering (CEMSE)
Numerical Porous Media Center (NumPor)
King Abdullah University of Science and Technology (KAUST)
http://numpor.kaust.edu.sa/

4700 King Abdullah University of Science and Technology
al-Khawarizmi Bldg (Bldg 1), Office # 4332
Thuwal 23955-6900, Kingdom of Saudi Arabia
http://www.kaust.edu.sa

Office Phone: +966 12 808-0459


Re: [petsc-users] TimeStepper norm problems. EMIL Please read this

2015-03-24 Thread Lisandro Dalcin
On 24 March 2015 at 06:15, Emil Constantinescu emcon...@mcs.anl.gov wrote:
 On 3/23/15 4:44 AM, Lisandro Dalcin wrote:

 Lisandro, that's a neat idea. If you are basically moving in the multistep
 realm for error estimation with one-step methods,

I'm wondering if no one ever did it. I could not find any reference.

 I think there are two
 considerations:

 1. The equations need to be continuous (no restarting, no events) or they
 would have to reset the error estimator somehow.


Indeed. My generalized-alpha code does not handle events at all. I'll
take a look into this. About restarting (not sure if we are talking
about the same thing), I have some code that is specific to
generalized-alpha (In the initial time step, I solve twice with
backward-Euler using dt/2 and dt) that let me properly initialize
generalized-alpha and also compute an initial error estimator.

 2. Backward differences are sensitive to time stepping changes; unlike
 one-step methods that can accommodate any step changes, B-D has strict
 limits for stability. While this does not play a crucial role (it is just
 the estimator);

Of course, though my B-D is order 2, and as you said, it is the just
for the estimator.

 if the time step varies wildly (can be easily constrained),
 it may create some problems with the estimator.


Indeed, the default clipping 0.1,10 is usually too wide, I usually do
-ts_adapt_basic_clip 0.5,2.0

Emil, is there any chance you can try my code with your problem? I
really need some feedback to push this to PETSc, otherwise

-- 
Lisandro Dalcin

Research Scientist
Computer, Electrical and Mathematical Sciences  Engineering (CEMSE)
Numerical Porous Media Center (NumPor)
King Abdullah University of Science and Technology (KAUST)
http://numpor.kaust.edu.sa/

4700 King Abdullah University of Science and Technology
al-Khawarizmi Bldg (Bldg 1), Office # 4332
Thuwal 23955-6900, Kingdom of Saudi Arabia
http://www.kaust.edu.sa

Office Phone: +966 12 808-0459


Re: [petsc-users] TimeStepper norm problems. EMIL Please read this

2015-03-24 Thread Lisandro Dalcin
On 24 March 2015 at 18:40, Emil Constantinescu emcon...@mcs.anl.gov wrote:
 Of course, though my B-D is order 2, and as you said, it is the just
 for the estimator.

 if the time step varies wildly (can be easily constrained),
 it may create some problems with the estimator.


 Indeed, the default clipping 0.1,10 is usually too wide, I usually do
 -ts_adapt_basic_clip 0.5,2.0


 The limits are a bit more strict for BDF2, but that's for the worst case and
 no filtering. So that looks good too.

I'm not sure I'm following you 100%. Please note I'm not using BDF2 in
the usual sense, what I do is to compute an A POSTERIORI error
estimator using a BD of order 2 to estimate LTE^{n+1} with the
solution vectors U^{n+1}, U^{n} and U^{n-1}, and this is AFTER
computing U^{n-1}.


-- 
Lisandro Dalcin

Research Scientist
Computer, Electrical and Mathematical Sciences  Engineering (CEMSE)
Numerical Porous Media Center (NumPor)
King Abdullah University of Science and Technology (KAUST)
http://numpor.kaust.edu.sa/

4700 King Abdullah University of Science and Technology
al-Khawarizmi Bldg (Bldg 1), Office # 4332
Thuwal 23955-6900, Kingdom of Saudi Arabia
http://www.kaust.edu.sa

Office Phone: +966 12 808-0459


Re: [petsc-users] TimeStepper norm problems. EMIL Please read this

2015-03-23 Thread Lisandro Dalcin
On 21 March 2015 at 17:32, Emil Constantinescu emcon...@mcs.anl.gov wrote:
 When -ts_theta_adapt is used, then it detects the instability as an error
 and reduces the step by a lot! wlte=1.24e+03 which means that the reduction
 should be severe but the controller tries 0.1*dt and that seems to pass but
 it jig-saws (take a look at the next attempted step), which means that it
 is likely unstable.

I think -ts_theta_adapt is seriously broken, I cannot make sense of
the way the error estimator is computed.

I have an alternative implementation of adaptivity in PetIGA in custom
Generalized-Alpha implementation. The adaptivity is not based on
embed-RK like for TSCN, but in storing a solution vector from the
previous step then estimate the LTE with backward differences. BTW,
this idea can be generalized to any second-order method at the cost of
tracking an extra solution vector and a bunch of Vec operations.

Emil, if you want to try my code to give it a try, it is very easy to
integrate, you just need to copy a file to your sources
(https://bitbucket.org/dalcinl/petiga/src/tip/src/tsalpha1.c) and add
a TSRegister() call after PetscInitialize(). Then you pass in the
command line -ts_type alpha1 -ts_alpha_adapt. You can also pass
-ts_alpha_radius 0.5 to add some high-frequency dissipation (the
default is 1.0, then the method is roughly equivalent to the midpoint
rule).

PS: I would like to put all this in PETSc, but before that, I need
help from you guys to make sure this method really makes sense. The
fact that this has worked very well for me in some problems is not
enough.

-- 
Lisandro Dalcin

Research Scientist
Computer, Electrical and Mathematical Sciences  Engineering (CEMSE)
Numerical Porous Media Center (NumPor)
King Abdullah University of Science and Technology (KAUST)
http://numpor.kaust.edu.sa/

4700 King Abdullah University of Science and Technology
al-Khawarizmi Bldg (Bldg 1), Office # 4332
Thuwal 23955-6900, Kingdom of Saudi Arabia
http://www.kaust.edu.sa

Office Phone: +966 12 808-0459


Re: [petsc-users] Memory leak in PetscRandom?

2015-02-26 Thread Lisandro Dalcin
On 26 February 2015 at 01:09, Jed Brown j...@jedbrown.org wrote:
 Barry Smith bsm...@mcs.anl.gov writes:

Thanks for the report. I ran the example with those options under
Linux with valgrind and found no memory leaks, I cannot run
valgrind on my Mac.

 When are you going to get a real operating system?  In the mean time,
 several people have reported that this works:

 http://ranf.tl/2014/11/28/valgrind-on-mac-os-x-10-10-yosemite/

When is valgrind going to use a real VCS? ;-)

-- 
Lisandro Dalcin

Research Scientist
Computer, Electrical and Mathematical Sciences  Engineering (CEMSE)
Numerical Porous Media Center (NumPor)
King Abdullah University of Science and Technology (KAUST)
http://numpor.kaust.edu.sa/

4700 King Abdullah University of Science and Technology
al-Khawarizmi Bldg (Bldg 1), Office # 4332
Thuwal 23955-6900, Kingdom of Saudi Arabia
http://www.kaust.edu.sa

Office Phone: +966 12 808-0459


Re: [petsc-users] Libraries destination directory is non-standard

2015-01-11 Thread Lisandro Dalcin
On 8 January 2015 at 23:32, Juan Luis Cano juanlu...@gmail.com wrote:
 Hello all,

 I am trying to build a conda package of PETSc 3.4 using its Python build
 system so anybody can install it in the Anaconda Python distribution. I need
 PETSc as a dependency for FEniCS.


Nice to know! I was thinking about doing it myself. You may find
interesting my mpich and openmpi packages (newer versions than the
available in Anaconda) to depend on https://binstar.org/mpi4py, I
still have to put the recipes in github repo.

 The problem is that the required shared libraries (also those corresponding
 to external packages such as UMFPACK) are installed in
 $PREFIX/lib/python2.7/site-packages/petsc/lib. The build process goes
 perfectly fine but when packaging this into a tarball these libs won't be
 found by any package as they are not in $PREFIX/lib. Here is the build
 script I am using:


Are you using pip to create/build the conda recipes? I really advice
against that. You should create a petsc package that is independent of
a Python runtime (of course, you still need Python to run configure),
this way you can install in the standard $PREFIX location. Why do you
say things are installed in site-packages? Who is installing there
FEniCS ?

Please note/remember that conda is not only for Python stuff, but for
any other Python-independent package that install under a $PREFIX
tree.

 https://github.com/Juanlu001/conda-recipes/blob/juanlu001/fenics/petsc/build.sh


* Please remember to add the md5 hash under the source section.
* python should not be a runtime dependency of the petsc package

 I tried to change the destination using `--prefix` in
 PETSC_CONFIGURE_OPTIONS and after `setup.py install`, but neither worked.
 Changing LD_LIBRARY_PATH seems dangerous to me and creating symbolic links
 three levels higher seems dirty. If the python version changes, or lib64
 is added somewhere, this will surely break... Does anybody have a suggestion
 on how to do this?


Knowing the full stack of conda packages you depend on would help to
figure out your issue.


-- 
Lisandro Dalcin

Research Scientist
Computer, Electrical and Mathematical Sciences  Engineering (CEMSE)
Numerical Porous Media Center (NumPor)
King Abdullah University of Science and Technology (KAUST)
http://numpor.kaust.edu.sa/

4700 King Abdullah University of Science and Technology
al-Khawarizmi Bldg (Bldg 1), Office # 4332
Thuwal 23955-6900, Kingdom of Saudi Arabia
http://www.kaust.edu.sa

Office Phone: +966 12 808-0459


Re: [petsc-users] Changing TSAdapt

2014-12-15 Thread Lisandro Dalcin
On 12 December 2014 at 01:44, Barry Smith bsm...@mcs.anl.gov wrote:

Miguel,

 Thanks for reporting this, you have found a bug in our code. When we 
 changed the adapt type we did not zero out the function pointers for the old 
 basic adaptor hence they were improperly called when the object was finally 
 destroyed at the end.

I've attached a patch. Once you apply this simply run

 make gnumake

in the PETSc root directory, recompile your code and run it again and it 
 should successfully end.

   Barry


Your patch also clears any user-specified adapt-ops-checkstage.
Perhaps you should revert your commit and merge my PR instead?

https://bitbucket.org/petsc/petsc/pull-request/228/fixes-for-ts-tsadapt-and-tsalpha/diff



-- 
Lisandro Dalcin

Research Scientist
Computer, Electrical and Mathematical Sciences  Engineering (CEMSE)
Numerical Porous Media Center (NumPor)
King Abdullah University of Science and Technology (KAUST)
http://numpor.kaust.edu.sa/

4700 King Abdullah University of Science and Technology
al-Khawarizmi Bldg (Bldg 1), Office # 4332
Thuwal 23955-6900, Kingdom of Saudi Arabia
http://www.kaust.edu.sa

Office Phone: +966 12 808-0459


Re: [petsc-users] Checkpointing / restart

2014-12-09 Thread Lisandro Dalcin
On 9 December 2014 at 18:24, Carles Bona carlesb...@gmail.com wrote:
 TS type is alpha


I'm working on a reimplementation of TSALPHA in PetIGA,
https://bitbucket.org/dalcinl/petiga/src/default/src/tsalpha1.c .
Eventually, I'll update the code in PETSc with this new
implementation. This new implementation features some attempts to
implement time-step adaptivity. Additionally, it implements a poor
man's attempt to estimate a initial derivative without requiring users
to setup and solve a problem involving the mass matrix.

We can experiment in this code with some extra APIs to let users
specify the initial derivative. What about TSAlphaSetSolution(ts,U,V)
and TSAlphaGetSolution(ts,U,V), where U and V are the initial
solution and derivative vectors?

-- 
Lisandro Dalcin

Research Scientist
Computer, Electrical and Mathematical Sciences  Engineering (CEMSE)
Numerical Porous Media Center (NumPor)
King Abdullah University of Science and Technology (KAUST)
http://numpor.kaust.edu.sa/

4700 King Abdullah University of Science and Technology
al-Khawarizmi Bldg (Bldg 1), Office # 4332
Thuwal 23955-6900, Kingdom of Saudi Arabia
http://www.kaust.edu.sa

Office Phone: +966 12 808-0459


Re: [petsc-users] PETSc/petsc4py error detection

2014-08-26 Thread Lisandro Dalcin
On 26 August 2014 13:24, Nicola Creati ncre...@inogs.it wrote:

 Thank you, but PetscPopSignalHandler has not been wrapped in petsc4py.


What about overriding the signal handler using Python's builtin signal module?


-- 
Lisandro Dalcin

Research Scientist
Computer, Electrical and Mathematical Sciences  Engineering (CEMSE)
Numerical Porous Media Center (NumPor)
King Abdullah University of Science and Technology (KAUST)
http://numpor.kaust.edu.sa/

4700 King Abdullah University of Science and Technology
al-Khawarizmi Bldg (Bldg 1), Office # 4332
Thuwal 23955-6900, Kingdom of Saudi Arabia
http://www.kaust.edu.sa

Office Phone: +966 12 808-0459


Re: [petsc-users] petsc4py getVecArray magic

2014-08-24 Thread Lisandro Dalcin
On 23 August 2014 05:15, Gianluca Meneghello gianm...@gmail.com wrote:
 Dear all,

 I just discovered that

 x = da.getVecArray(X)

 does not return a numpy array but rather a _DMDA_Vec_array, which has (I
 guess) the great advantage of allowing for a starting index which is not
 zero. So x[87] makes sense even if the actual length of x is only 50.

 What I would like to do is to split the _DMDA_Vec_array x into multiple
 variables, e.g.

 x1,x2,x3 = x[:,0:12] , x[:,12:16], x[:,16:28]

 which, written in this way, returns x1,x2,x3 as numpy view() arrays starting
 from zero and of length equal to the length of the global X divided by the
 number of processors. The arrays x1,x2,x3 thus becomes unaccessible when
 using global indices. Using local indices does not look like a good idea
 because of the presence of ghost entries.

 I was wondering if there is any good way to address this situation besides
 not splitting x in x1,x2,x3.


Sorry, but I'm not following you.


-- 
Lisandro Dalcin

Research Scientist
Computer, Electrical and Mathematical Sciences  Engineering (CEMSE)
Numerical Porous Media Center (NumPor)
King Abdullah University of Science and Technology (KAUST)
http://numpor.kaust.edu.sa/

4700 King Abdullah University of Science and Technology
al-Khawarizmi Bldg (Bldg 1), Office # 4332
Thuwal 23955-6900, Kingdom of Saudi Arabia
http://www.kaust.edu.sa

Office Phone: +966 12 808-0459


Re: [petsc-users] MPI-IO

2014-07-22 Thread Lisandro Dalcin
On 21 July 2014 13:50, Jed Brown j...@jedbrown.org wrote:

 Must have been prior to MPI-2.1.  Shall we test for MPI_Pack_external
 and use MPI-IO by default if it's available?

Be prepared to file tickets for MPICH and OpenMPI.
Pack/unpack_external with MPI_DOUBLE never worked for me:
https://bitbucket.org/mpi4py/mpi4py/src/master/test/test_pack.py#cl-105


-- 
Lisandro Dalcin
---
CIMEC (UNL/CONICET)
Predio CONICET-Santa Fe
Colectora RN 168 Km 472, Paraje El Pozo
3000 Santa Fe, Argentina
Tel: +54-342-4511594 (ext 1016)
Tel/Fax: +54-342-4511169


Re: [petsc-users] Vector Operations

2013-12-26 Thread Lisandro Dalcin
On 25 December 2013 11:43, 吕超 luc...@mail.iggcas.ac.cn wrote:
 So, I want to kown how can I  get c = a.*b in Petsc. Can you help me?


http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Vec/VecPointwiseMult.html

-- 
Lisandro Dalcin
---
CIMEC (UNL/CONICET)
Predio CONICET-Santa Fe
Colectora RN 168 Km 472, Paraje El Pozo
3000 Santa Fe, Argentina
Tel: +54-342-4511594 (ext 1016)
Tel/Fax: +54-342-4511169


Re: [petsc-users] MatMPIBAIJSetPreallocationCSR i and j indices

2013-09-29 Thread Lisandro Dalcin
On 29 September 2013 17:31, Barry Smith bsm...@mcs.anl.gov wrote:

I forgot that MatSetValuesBlocked() takes row oriented input even within 
 the blocks. So if we simply call 
 MatSetOption(mat,MAT_ROW_ORIENTED,PETSC_FALSE) on the matrix after it is 
 created then the storage should work as I originally planned.

   Barry

 I found this discussion values(bs,nblockcols,bs,nblockrows)  very confusing. 
 MatMPIBAIJSetPreallocationCSR() doesn't take a square block of values so I 
 didn't see what you are talking about. You must be talking about each 
 internal call to 
 MatSetValuesBlocked_MPIBAIJ(B,1,row,ncols,icols,svals,INSERT_VALUES) where 
 it puts in one block row of values with each call.


Yes, that's what I was talking about. And yes, if you have
column-oriented block values [e.g, in fotran you use an array
VALUES(bs,bs,nnz) ] and use MatSetOption(), you are fine and
MatMPIBAIJSetPreallocationCSR() will do the right thing. Otherwise
(e.g, you are in C and use VALUES[nnz][bs][bs]), feeding
MatMPIBAIJSetPreallocationCSR() with these values will not work. I
guess we have to live with that.


-- 
Lisandro Dalcin
---
CIMEC (UNL/CONICET)
Predio CONICET-Santa Fe
Colectora RN 168 Km 472, Paraje El Pozo
3000 Santa Fe, Argentina
Tel: +54-342-4511594 (ext 1016)
Tel/Fax: +54-342-4511169


Re: [petsc-users] MatMPIBAIJSetPreallocationCSR i and j indices

2013-09-27 Thread Lisandro Dalcin
On 27 September 2013 21:23, Barry Smith bsm...@mcs.anl.gov wrote:
  Just use MatCreateMPIBAIJWithArrays() saves thinking and work.

BTW, Have you ever thought about how hard to use is this API if you
want to provide the matrix values? How are users supposed to fill the
values array for a matrix with bs1 ? This API is certainly not
Fortran-friendly (and I think it not even C-friendly). Mateo had the
CSR indices as well as the block values, but I recommended to use the
CSR indices to perform proper matrix preallocation, then loop over
cells and call MatSetValuesBlocked(). Otherwise, he would need to
create an new intermediate array and fill it the right way for
MatCreateMPIBAIJWithArrays() to work as expected.

The root of the problem is how MatSetValuesBlocked() expects the array
of values to be layout in memory. While I understand that the current
convention make it compatible with MatSetValues(), I'm always
wondering if a layout like values[nblockrows][nblockcols][bs][bs]
would not be much more convenient for users.



-- 
Lisandro Dalcin
---
CIMEC (UNL/CONICET)
Predio CONICET-Santa Fe
Colectora RN 168 Km 472, Paraje El Pozo
3000 Santa Fe, Argentina
Tel: +54-342-4511594 (ext 1016)
Tel/Fax: +54-342-4511169


[petsc-users] petsc 3.4:

2013-05-15 Thread Lisandro Dalcin
On 15 May 2013 16:54, Satish Balay balay at mcs.anl.gov wrote:
 On Wed, 15 May 2013, Matteo Parsani wrote:

 Hello Satish,
 Lisandro has found that in the make file we have

 include ././${PETSC_ARCH}/conf/petscvariables

 we try to support too many ways of installing petsc [with prefix/
 without prefix, with a defaut PETSC_ARCH, without PETSC_DIR set etc..]
 And we need this line for users who forget to set PETSC_ARCH in a
 non-prefix build. [to pick up a defaut PETSC_ARCH]



Satish, I think all what is needed is to fix a little the test_build
target. Basically, test if the file under PETSC_ARCH exists, otherwise
use PETSC_DIR/include/petscconf.h . Perhpas I'm missing something, but
this should be near to work.



--
Lisandro Dalcin
---
CIMEC (INTEC/CONICET-UNL)
Predio CONICET-Santa Fe
Colectora RN 168 Km 472, Paraje El Pozo
3000 Santa Fe, Argentina
Tel: +54-342-4511594 (ext 1011)
Tel/Fax: +54-342-4511169


[petsc-users] petsc 3.4:

2013-05-15 Thread Lisandro Dalcin
On 15 May 2013 17:12, Satish Balay balay at mcs.anl.gov wrote:
 On Wed, 15 May 2013, Lisandro Dalcin wrote:

 On 15 May 2013 16:54, Satish Balay balay at mcs.anl.gov wrote:
  On Wed, 15 May 2013, Matteo Parsani wrote:
 
  Hello Satish,
  Lisandro has found that in the make file we have
 
  include ././${PETSC_ARCH}/conf/petscvariables
 
  we try to support too many ways of installing petsc [with prefix/
  without prefix, with a defaut PETSC_ARCH, without PETSC_DIR set etc..]
  And we need this line for users who forget to set PETSC_ARCH in a
  non-prefix build. [to pick up a defaut PETSC_ARCH]
 


 Satish, I think all what is needed is to fix a little the test_build
 target. Basically, test if the file under PETSC_ARCH exists, otherwise
 use PETSC_DIR/include/petscconf.h . Perhpas I'm missing something, but
 this should be near to work.

 We've avoided using gnu make extensions - so there is no usage of 'if'
 directives in petsc makefiles. [and we've used include directive as
 alternative]


I was not clear enough. I'm talking about a shell if, something like

if [-f ${PETSC_DIR}/${PETSC_ARCH}/include/petscconf.h]; then \
   egrep ...
else
   egrep 
fi


--
Lisandro Dalcin
---
CIMEC (INTEC/CONICET-UNL)
Predio CONICET-Santa Fe
Colectora RN 168 Km 472, Paraje El Pozo
3000 Santa Fe, Argentina
Tel: +54-342-4511594 (ext 1011)
Tel/Fax: +54-342-4511169


[petsc-users] Is there any counterpart of TSSetMatrices in petsc4py?

2012-08-15 Thread Lisandro Dalcin
On 15 August 2012 14:13, Xin Zhao sean.null at gmail.com wrote:
 Dear all,

 What I need is just to solve a very easy linear ODE.

 u_t = Au

 A is sparse and constant.

 But I didn't find a counterpart of TSSetMatrices in petsc4py for this.
 It is strange because petsc4py can handle much harder problems.


TSSetMatrices() is no longer supported in petsc-3.3. You have to
rewrite your problem in implicit nonlinear form, and ask TS to solve a
linear problem. You still have to code your IFunction and IJacobian in
Python, use the PETSc C implementations TSComputeIFunctionLinear() and
TSComputeIJacobianConstant() and translate to Python.


-- 
Lisandro Dalcin
---
CIMEC (INTEC/CONICET-UNL)
Predio CONICET-Santa Fe
Colectora RN 168 Km 472, Paraje El Pozo
3000 Santa Fe, Argentina
Tel: +54-342-4511594 (ext 1011)
Tel/Fax: +54-342-4511169


[petsc-users] Getting 0 for both start and end with MatGetOwnershipRange()

2012-05-09 Thread Lisandro Dalcin
On 9 May 2012 20:01, Daniel Christian Weflen Daniel.Weflen at colorado.edu 
wrote:
 I have the following code:

 ??? Mat H;
 ??? PetscInt mat_size=x_npts*y_npts;
 ??? ierr=MatCreate(PETSC_COMM_WORLD,H);CHKERRQ(ierr);

 ierr=MatSetSizes(H,PETSC_DECIDE,PETSC_DECIDE,mat_size,mat_size);CHKERRQ(ierr);
 ??? ierr=MatSetType(H,MATSEQMAIJ);CHKERRQ(ierr);
 ??? PetscInt* nonzeros;
 ??? PetscInt start_row,end_row,row_index;
 ??? MatGetOwnershipRange(H,start_row,end_row);

 When I check the values of start_row and end_row right after calling
 MatGetOwnershipRange, both are 0. What am I doing wrong here? Am I calling
 MatCreate(), MatSetSizes() and MatSetType() in the wrong order?

Call MatSetUp() after MatSetType().


-- 
Lisandro Dalcin
---
CIMEC (INTEC/CONICET-UNL)
Predio CONICET-Santa Fe
Colectora RN 168 Km 472, Paraje El Pozo
3000 Santa Fe, Argentina
Tel: +54-342-4511594 (ext 1011)
Tel/Fax: +54-342-4511169


[petsc-users] petsc4py and numpy

2012-05-08 Thread Lisandro Dalcin
On 8 May 2012 12:04, Matthew Knepley knepley at gmail.com wrote:
 On Tue, May 8, 2012 at 11:00 AM, Christian Staudt
 christian.staudt at ira.uka.de wrote:

 Matthew Knepley?knepley at gmail.com?wrote:

 a) Is there an easy and efficient way to convert a numpy.ndarray to a

 PETSc.Vec or PETSc.Mat (and vice-versa)? ?(A PETSc.Vec.getArray() returns
 a

 numpy.ndarray, though I have found no such method for PETSc.Mat yet. )



 If you are using dense matrices (MATDENSE) then MatGetArray() works.


 I don't think this is implemented in petsc4py: At least, there's no method
 PETSc.Mat.getArray()


 If you ask Lisandro for it, he will put it in.


Matt, I've never put that in just because I'm not sure how to make it
general. PETSc have many different matrix formats, and MatGetArray()
call is like a hack, where you can get very different things (the 1D
values of an AIJ matrix, the 2D values of a dense matrix in Fortran
order). What about MPIAIJ, what should I return to users? Can you
suggest some approach or API to implement this (I mean, how to get
arrays for the various matrix types)?

-- 
Lisandro Dalcin
---
CIMEC (INTEC/CONICET-UNL)
Predio CONICET-Santa Fe
Colectora RN 168 Km 472, Paraje El Pozo
3000 Santa Fe, Argentina
Tel: +54-342-4511594 (ext 1011)
Tel/Fax: +54-342-4511169


[petsc-users] petsc4py and numpy

2012-05-08 Thread Lisandro Dalcin
On 8 May 2012 12:25, Matthew Knepley knepley at gmail.com wrote:
 On Tue, May 8, 2012 at 11:15 AM, Lisandro Dalcin dalcinl at gmail.com wrote:

 On 8 May 2012 12:04, Matthew Knepley knepley at gmail.com wrote:
  On Tue, May 8, 2012 at 11:00 AM, Christian Staudt
  christian.staudt at ira.uka.de wrote:
 
  Matthew Knepley?knepley at gmail.com?wrote:
 
  a) Is there an easy and efficient way to convert a numpy.ndarray to a
 
  PETSc.Vec or PETSc.Mat (and vice-versa)? ?(A PETSc.Vec.getArray()
  returns
  a
 
  numpy.ndarray, though I have found no such method for PETSc.Mat yet. )
 
 
 
  If you are using dense matrices (MATDENSE) then MatGetArray() works.
 
 
  I don't think this is implemented in petsc4py: At least, there's no
  method
  PETSc.Mat.getArray()
 
 
  If you ask Lisandro for it, he will put it in.
 

 Matt, I've never put that in just because I'm not sure how to make it
 general. PETSc have many different matrix formats, and MatGetArray()
 call is like a hack, where you can get very different things (the 1D
 values of an AIJ matrix, the 2D values of a dense matrix in Fortran
 order). What about MPIAIJ, what should I return to users? Can you
 suggest some approach or API to implement this (I mean, how to get
 arrays for the various matrix types)?


 I think the whole point is that you should not do it. However, that has
 never stopped us from giving users
 what they ask for. We are not int he business of reforming computing
 practices, which is why we have
 MatGetArray(). I would make a wrapper that failed unless the type was
 MATDENSE.


Yea, that's really easy.

However, Is there anything better for AIJ matrices? Right now, we have
mat.getValuesCSR()... it works by calling MatGetRow() row by row twice
(first pass to count and allocate, next to get values and copy to
NumPy arrays).


-- 
Lisandro Dalcin
---
CIMEC (INTEC/CONICET-UNL)
Predio CONICET-Santa Fe
Colectora RN 168 Km 472, Paraje El Pozo
3000 Santa Fe, Argentina
Tel: +54-342-4511594 (ext 1011)
Tel/Fax: +54-342-4511169


[petsc-users] Petsc CMake conf

2012-04-14 Thread Lisandro Dalcin
On 14 April 2012 01:19, Jed Brown jedbrown at mcs.anl.gov wrote:


 On Fri, Apr 13, 2012 at 17:17, Mohammad Mirzadeh mirzadeh at gmail.com 
 wrote:

 Hi,

 Are there any cmake conf file that I can include in my CMakeLists.txt that
 would automatically detect include/lib directory for my petsc installation?
 Right now I'm doing this manually and was wondering if its possibel to
 automate it.


 https://github.com/jedbrown/cmake-modules/blob/master/FindPETSc.cmake


Jed, have you considered adding this to petsc-dev, removing all the
checks for old PETSc releases? Or perhaps better, how about  dumping
an autogenerated PETSc.cmake in
$PETSC_DIR/$PETSC_ARCH/conf/PETSc.cmake, or perhaps other name,
defining all the stuff end users need. For example, I'm using the code
below in a codebase I'm  working on (note that at the end I'm loading
PETScConfig.cmake), but this is not robust, it does not handle builds
with separate libraries. If PETSc provides all these variables, then I
could simply  include (${PETSC_DIR}/${PETSC_ARCH}/conf/PETSc.cmake),
and you can make your FindPETSc.cmake much simpler (at least for the
upcoming version of PETSc)


  find_path (PETSC_DIR include/petsc.h HINTS ENV PETSC_DIR DOC PETSc
Directory)
  set (PETSC_ARCH $ENV{PETSC_ARCH})
  find_path (PETSC_INCLUDE_DIR  petsc.h HINTS ${PETSC_DIR}
PATH_SUFFIXES include NO_DEFAULT_PATH)
  find_path (PETSC_INCLUDE_CONF petscconf.h HINTS ${PETSC_DIR}
PATH_SUFFIXES ${PETSC_ARCH}/include NO_DEFAULT_PATH)
  mark_as_advanced (PETSC_INCLUDE_DIR PETSC_INCLUDE_CONF)
  set (PETSC_INCLUDES ${PETSC_INCLUDE_CONF} ${PETSC_INCLUDE_DIR} CACHE
PATH PETSc include paths FORCE)
  set (PETSC_DEFINITIONS -D__INSDIR__= CACHE STRING PETSc
preprocesor definitions FORCE)
  find_path (PETSC_LIB_DIR NAMES  HINTS ${PETSC_DIR} PATH_SUFFIXES
${PETSC_ARCH}/lib NO_DEFAULT_PATH)
  find_library (PETSC_LIBRARY NAMES petsc HINTS ${PETSC_LIB_DIR}
NO_DEFAULT_PATH)
  set (PETSC_LIBRARIES ${PETSC_LIBRARY} CACHE FILEPATH PETSc library FORCE)
  include (${PETSC_DIR}/${PETSC_ARCH}/conf/PETScConfig.cmake)


-- 
Lisandro Dalcin
---
CIMEC (INTEC/CONICET-UNL)
Predio CONICET-Santa Fe
Colectora RN 168 Km 472, Paraje El Pozo
3000 Santa Fe, Argentina
Tel: +54-342-4511594 (ext 1011)
Tel/Fax: +54-342-4511169


[petsc-users] petsc4py KSP Question

2012-03-07 Thread Lisandro Dalcin
On 7 March 2012 02:22, Gaetan Kenway kenway at utias.utoronto.ca wrote:
 Hello

 I'm in the process of using petsc4py to solve a large multidisciplinary,
 non-linear system and its adjoint. I have the non-linear solution with
 snes() working correctly and I'm now doing the linear solution.

 For the non-linear solve, I create the snes and set my user context as
 follows:
 ? ?# Create SNES Object
 ? ASContext = ASNKSolver(self) # Context to hold data
 ? snes = PETSc.SNES().createPython(ASContext, comm=self.gcomm)
 ? snes.setFunction(ASContext.formFunction, resVec)

 For the linear part, I'm a little confused. I currently have

 ? ? # Create Python Context and KSP Object
 ? ? KSPContext = AdjointKSPSolver(self) # Context to hold data
 ? ? ksp = PETSc.KSP().createPython(KSPContext, comm=self.gcomm)

 What I'm not sure of is how to specify the operator for the KSP solver in
 the KSPContext object. Is this possible? Is there something like def apply()
 you must do?

 The petsc4py poisson2d.py example, first creates a Mat, then sets a context
 for that, and then uses ksp.setOperators() to use that matrix. Is this the
 only way to do it? If, so, what is the use of
 the??PETSc.KSP().createPython() command?


Unless you want to implement a custom KSP/SNES in Python, you should
just use KSP/SNES.create(). KSP/SNES.createPython() is a rather
advanced feature, it is not obvious at all how to use them, there are
no demos nor documentation about it (however, there is some code in
test/test_{ksp|snes}_py.py ).

Please take a look at the demo/bratu2d and demo/bratu3d about how to
setup a SNES solver. For KSP, take a look at
demo/kspsolve/petsc-mat.py and demo/kspsolve/petsc-ksp.py



-- 
Lisandro Dalcin
---
CIMEC (INTEC/CONICET-UNL)
Predio CONICET-Santa Fe
Colectora RN 168 Km 472, Paraje El Pozo
3000 Santa Fe, Argentina
Tel: +54-342-4511594 (ext 1011)
Tel/Fax: +54-342-4511169


[petsc-users] invalid object classid 0 between PETSc and SLEPc

2011-12-28 Thread Lisandro Dalcin
On 28 December 2011 01:16, Geoffrey Irving irving at naml.us wrote:
 Hello,

 I've installed PETSc 3.2 and SLEPc 3.2 on Mac OS 10.7 through MacPorts
 2.0.3. ?PETSc works fine by itself, but SLEPc fails on the initial
 EPSCreate call. ?The error follows.

 This is presumably a MacPorts misconfiguration issue, but thought I'd
 check here first to see if anyone has seen this before in a similar
 situation. ?Only static libraries exist, so as far as I know the
 suggested explanation doesn't apply. ?The conf directory is at
 http://naml.us/~irving/random/petsc-conf.tar.gz if it helps.


Just in case, please run:

$ otool -L ./your_executable

Do you have any other PETSc/SLEPc build lying around?

-- 
Lisandro Dalcin
---
CIMEC (INTEC/CONICET-UNL)
Predio CONICET-Santa Fe
Colectora RN 168 Km 472, Paraje El Pozo
3000 Santa Fe, Argentina
Tel: +54-342-4511594 (ext 1011)
Tel/Fax: +54-342-4511169


[petsc-users] petsc4py: PetSc is no longer initialised when calling C Function from Cython

2011-10-21 Thread Lisandro Dalcin
On 20 October 2011 19:34, Jeff Wiens jeffrey.k.wiens at gmail.com wrote:

 Thanks for your help. It saved me a LOT of time.
 Jeff


BTW, have you ever used SWIG? If the functions you need to wrap are
simple (let say, any PetscObject subtype and scalar paramenters) you
can get your wrappers with less code to write on your side.


-- 
Lisandro Dalcin
---
CIMEC (INTEC/CONICET-UNL)
Predio CONICET-Santa Fe
Colectora RN 168 Km 472, Paraje El Pozo
3000 Santa Fe, Argentina
Tel: +54-342-4511594 (ext 1011)
Tel/Fax: +54-342-4511169


[petsc-users] petsc4py: PetSc is no longer initialised when calling C Function from Cython

2011-10-20 Thread Lisandro Dalcin
On 20 October 2011 17:05, Jeff Wiens jeffrey.k.wiens at gmail.com wrote:

 The problem seems to be that PetSc is no longer initialised when you
 call a C function. If I initialise PETSc from C, the program works.
 However, if I don't call init() from C, the program won't use the
 python PetSc.

Mmm... I bet you built PETSc with static libraries. There is not
(easyportable) way to mix petsc4py and other C codes when using
static libraries. Please build PETSc with shared libraries.


 Eventually, I would like
 to Initialise PETSc and construct my PETSc vectors in python and then
 have my C code perform operations on them. Again, this yet another
 test problem.


Have you seen the demo/wrap-cython example in petsc4py sources?


-- 
Lisandro Dalcin
---
CIMEC (INTEC/CONICET-UNL)
Predio CONICET-Santa Fe
Colectora RN 168 Km 472, Paraje El Pozo
3000 Santa Fe, Argentina
Tel: +54-342-4511594 (ext 1011)
Tel/Fax: +54-342-4511169


[petsc-users] petsc4py: PetSc is no longer initialised when calling C Function from Cython

2011-10-20 Thread Lisandro Dalcin
On 20 October 2011 19:34, Jeff Wiens jeffrey.k.wiens at gmail.com wrote:
 Lisandro, the shared library was the problem. I seem to have
 everything working by including the --with-shared and --with-dynamic
 options in my PetSc 3.1 configuration.


The --with-dynamic stuff should not be required, and I think it's
better if you remove it.

 Note, I had to make significant modifications to your setup.py file
 for your cython demo to work (on my installation).


Expected

 I have attached my
 setup.py file, which is based on code I found in PyClaw. The main
 differences is overloading:
 ? ? ? ? build_src.build_src.generate_a_pyrex_source = generate_a_cython_source
 and including the following directories:
 ? ?INCLUDE_DIRS +=[/PATH/To/openmpi/1.4.3/include]
 ? ?INCLUDE_DIRS += [/PATH/To/hdf5/1.8.5-patch1/include]

 Thanks for your help. It saved me a LOT of time.
 Jeff


I'll take a look at your setup.py, thanks!


-- 
Lisandro Dalcin
---
CIMEC (INTEC/CONICET-UNL)
Predio CONICET-Santa Fe
Colectora RN 168 Km 472, Paraje El Pozo
3000 Santa Fe, Argentina
Tel: +54-342-4511594 (ext 1011)
Tel/Fax: +54-342-4511169


[petsc-users] Naturally ordered to DA vector copy

2011-07-20 Thread Lisandro Dalcin
On 20 July 2011 01:34, khalid ashraf khalid_eee at yahoo.com wrote:
 Hello,
 In my application, I need to do both Finite element(FEM) and Finite
 Difference(FD) computations.
 I am doing the FD using the DA and FEM using natural numbering(x changes
 fastest, then y then z).
 My question is what would be the fastest way to copy the FEM vectors to the
 FD vectors. Looks like
 DAGetAO may provide a mapping between DA vector layout and natural layout.
 Since I am not very familiar with programming using AO, just wanted to ask
 if there is a more efficient/easier way to copy from the naturally ordered
 vector to the DA vector before delving into AO. If I have to use AO, could
 you please redirect me to some related example.
 Thanks in advance.
 Khalid

DMDAGlobalToNaturalBegin/End and DMDANaturalToGlobalBegin/End

-- 
Lisandro Dalcin
---
CIMEC (INTEC/CONICET-UNL)
Predio CONICET-Santa Fe
Colectora RN 168 Km 472, Paraje El Pozo
3000 Santa Fe, Argentina
Tel: +54-342-4511594 (ext 1011)
Tel/Fax: +54-342-4511169


[petsc-users] PETSc recommended visualization packages

2011-07-14 Thread Lisandro Dalcin
On 14 July 2011 16:22, Ethan Coon ecoon at lanl.gov wrote:
 Attached patch takes Lisandro's matio.py and puts it in a
 PetscBinaryRead.py that reads a binary with (potentially) multiple Petsc
 objects, deciphers the contents of the file from the header, and returns
 numpy objects. ?It's basically the same as the PetscBinaryRead.m

 I've tested this on Vecs, but not Mats or ISs, and I haven't thoroughly
 tested degenerate cases. ?Will do some testing, but I wanted Lisandro to
 take a look at it as well...


1) readMatDense() is wrong. The I,J,V arrays are the CSR structure,
not the COO (coordinate) format, Then you cannot simply:

   for i,j,v in zip(I,J,V):
mat[i,j] = v

I think you have to:

for row, rstart in enumerate(I):
rend = I[row+1]
mat[row, J[rstart:rend]] = V[rstart:rend]


2) It would be nice to provie a 'scipy' Mat format returning a
'scipy.sparse.csr_matrix' instance. A possible implementation would
be:

def readMatSciPy(fh):
(M, N), (I, J, V) = readMatSparse(fh)
from scipy.sparse import csr_matrix
return csr_matrix((V, J, I), shape=(M, N))



-- 
Lisandro Dalcin
---
CIMEC (INTEC/CONICET-UNL)
Predio CONICET-Santa Fe
Colectora RN 168 Km 472, Paraje El Pozo
3000 Santa Fe, Argentina
Tel: +54-342-4511594 (ext 1011)
Tel/Fax: +54-342-4511169


[petsc-users] PETSc recommended visualization packages

2011-07-14 Thread Lisandro Dalcin
On 14 July 2011 18:33, Ethan Coon ecoon at lanl.gov wrote:
 On Thu, 2011-07-14 at 17:47 -0300, Lisandro Dalcin wrote:


 1) readMatDense() is wrong. The I,J,V arrays are the CSR structure,
 not the COO (coordinate) format, Then you cannot simply:


 Ok, thanks, I wasn't sure of which format you were using,

Well, PETSc's AIJ matrices are CSR in memory (but not in disk, binary
files store row counts instead of  row pointers)


 2) It would be nice to provie a 'scipy' Mat format returning a
 'scipy.sparse.csr_matrix' instance. A possible implementation would
 be:


 Agreed, I considered doing that as the default for sparse, but wasn't
 sure what others would want. ?Providing both makes sense.


Note however that this would require a SciPy installation.

PS: Sorry for not being more helpful and contributing actual code, but
I'm really busy.


-- 
Lisandro Dalcin
---
CIMEC (INTEC/CONICET-UNL)
Predio CONICET-Santa Fe
Colectora RN 168 Km 472, Paraje El Pozo
3000 Santa Fe, Argentina
Tel: +54-342-4511594 (ext 1011)
Tel/Fax: +54-342-4511169


[petsc-users] PETSc recommended visualization packages

2011-07-06 Thread Lisandro Dalcin
On 6 July 2011 11:41, Ethan Coon ecoon at lanl.gov wrote:
 Both Paraview and Visit generally make ugly axes/colorbars/keys/etc. ?In my 
 opinion they both look fine for presentations and my own viewing, but are not 
 really acceptable for publication-quality. ?For things run on a 
 reasonably-sized problem, matplotlib (using either h5py to read hdf5, pyvtk 
 to read vtk, or petsc4py to read PETSc's binary format) makes publishable 
 plots.


Using petsc4py for reading PETSc binary formats is overkill. All this
can be done with just NumPy. Should PETSc ship some Python code to
read IS's, Vec's, Mat's in binary format? Can any of you suggest an
appropriate location in the source tree?




-- 
Lisandro Dalcin
---
CIMEC (INTEC/CONICET-UNL)
Predio CONICET-Santa Fe
Colectora RN 168 Km 472, Paraje El Pozo
3000 Santa Fe, Argentina
Tel: +54-342-4511594 (ext 1011)
Tel/Fax: +54-342-4511169


[petsc-users] PETSc recommended visualization packages

2011-07-05 Thread Lisandro Dalcin
On 5 July 2011 12:53, Brad Aagaard baagaard at usgs.gov wrote:
 Hi,

 I agree with Blaise's assessment of ParaView. It is good for quick and
 dirty visualization, but its scripting interface is very limited.

 MayaVi provides a GUI along with high-level API access to VTK data
 structures through Python. One can update/modify the data values through
 the API without ever writing VTK compatible files. I read my data from
 HDF5 files using PyTables and create VTK data objects using the API. A
 similar mechanism could be used for data stored using SQLite. I think
 some Linux and Mac packaging systems include MayaVi; I usually build
 from source.


I recently added support (requires petsc4py, of course) for
-ksp|snes|ts_monitor_python filename.py:function|class. Perhaps you
can find this useful for live monitoring of simulations.


-- 
Lisandro Dalcin
---
CIMEC (INTEC/CONICET-UNL)
Predio CONICET-Santa Fe
Colectora RN 168 Km 472, Paraje El Pozo
3000 Santa Fe, Argentina
Tel: +54-342-4511594 (ext 1011)
Tel/Fax: +54-342-4511169


[petsc-users] Question for time dependent problems

2011-07-05 Thread Lisandro Dalcin
On 5 July 2011 17:35, Barry Smith bsmith at mcs.anl.gov wrote:

 Beginning of a time step (1st implicit stage of RK scheme):
 1. Use MatZeroEntries on the preconditioner

 ?Do not call MatZeroEntries on a freshly created matrix (that destroys the 
 preallocation pattern) so skip the MatZeroEntries the first time.


What do you mean? The OP said that the nonzero structure was set at
problem set up.

-- 
Lisandro Dalcin
---
CIMEC (INTEC/CONICET-UNL)
Predio CONICET-Santa Fe
Colectora RN 168 Km 472, Paraje El Pozo
3000 Santa Fe, Argentina
Tel: +54-342-4511594 (ext 1011)
Tel/Fax: +54-342-4511169


[petsc-users] petsc4py DA with dof1 confusion

2011-07-04 Thread Lisandro Dalcin
On 4 July 2011 14:25, Matthew Emmett memmett at unc.edu wrote:
 Hi everyone,

 I am trying to use petsc4py in a fairly straight forward finite volume
 shallow-water code (the PDE has three unknowns). I have a structured
 grid and am using a DA to communicate ghost cells. ?When I use dof=1,
 everything works as I would expect. ?However, when I try to use dof=3,
 I am confused by the results. ?I have probably mis-understood
 something -- any suggestions would be appreciated.

 For example (see attached script), first I create the DA ('dof' and
 'width' are defined):

 ?da = PETSc.DA().create(dim=2, sizes=[8, 8], dof=dof,
 ? ? ? ? ? ? ? ? ? ? ? ? boundary_type='periodic',
 ? ? ? ? ? ? ? ? ? ? ? ? stencil_width=width,
 ? ? ? ? ? ? ? ? ? ? ? ? stencil_type='box')

 Next, I create global and local vectors:

 ?gvec = da.createGlobalVec()
 ?lvec = da.createLocalVector()

 Then, put something into the local vectors, and scatter:

 ?with da.getVecArray(lvec) as va:
 ? ?if dof == 1:
 ? ? ?va.array[width:-width,width:-width] = MPI.COMM_WORLD.rank+1
 ? ?else:
 ? ? ?va.array[width:-width,width:-width,0] = MPI.COMM_WORLD.rank+1
 ?da.localToGlobal(lvec, gvec)

 Finally, get the local vector again and print:

 ?da.globalToLocal(gvec, lvec)
 ?with da.getVecArray(lvec) as va:
 ? ?if dof == 1:
 ? ? ?print va.array[:,:]
 ? ?else:
 ? ? ?print va.array[:,:,0]

 With 'dof=1' I get the following on the second processor:

 [[ 1. ?2. ?2. ?2. ?2. ?1.]
 ?[ 1. ?2. ?2. ?2. ?2. ?1.]
 ?[ 1. ?2. ?2. ?2. ?2. ?1.]
 ?[ 1. ?2. ?2. ?2. ?2. ?1.]
 ?[ 1. ?2. ?2. ?2. ?2. ?1.]
 ?[ 1. ?2. ?2. ?2. ?2. ?1.]
 ?[ 1. ?2. ?2. ?2. ?2. ?1.]
 ?[ 1. ?2. ?2. ?2. ?2. ?1.]
 ?[ 1. ?2. ?2. ?2. ?2. ?1.]
 ?[ 1. ?2. ?2. ?2. ?2. ?1.]]

 With 'dof=3' I get the following on the second processor:

 [[ 0. ?0. ?0. ?0. ?0. ?0.]
 ?[ 0. ?0. ?0. ?0. ?2. ?0.]
 ?[ 0. ?0. ?0. ?0. ?2. ?0.]
 ?[ 0. ?0. ?0. ?2. ?2. ?0.]
 ?[ 0. ?0. ?0. ?2. ?2. ?0.]
 ?[ 0. ?0. ?0. ?2. ?2. ?0.]
 ?[ 0. ?0. ?0. ?2. ?2. ?0.]
 ?[ 0. ?0. ?0. ?2. ?2. ?2.]
 ?[ 0. ?0. ?0. ?2. ?2. ?2.]
 ?[ 0. ?0. ?0. ?0. ?0. ?2.]]

 I was expecting to see the same as the 'dof=1' case, since I am
 setting (and printing) va.array[:,:,0].


 Again, any suggestions would be appreciated. ?Thanks,
 Matt

 PS, thanks to the developers of PETSc and petsc4py for all your hard
 work. ?I had the pleasure of meeting Jed and Matt at KAUST recently.


Mmm, it seems that this functionality is broken for the case of
boundary_type='periodic'. I'll take a look.


-- 
Lisandro Dalcin
---
CIMEC (INTEC/CONICET-UNL)
Predio CONICET-Santa Fe
Colectora RN 168 Km 472, Paraje El Pozo
3000 Santa Fe, Argentina
Tel: +54-342-4511594 (ext 1011)
Tel/Fax: +54-342-4511169


[petsc-users] petsc4py DA with dof1 confusion

2011-07-04 Thread Lisandro Dalcin
On 4 July 2011 15:17, Lisandro Dalcin dalcinl at gmail.com wrote:
 On 4 July 2011 14:25, Matthew Emmett memmett at unc.edu wrote:
 Hi everyone,

 I am trying to use petsc4py in a fairly straight forward finite volume
 shallow-water code (the PDE has three unknowns). I have a structured
 grid and am using a DA to communicate ghost cells. ?When I use dof=1,
 everything works as I would expect. ?However, when I try to use dof=3,
 I am confused by the results. ?I have probably mis-understood
 something -- any suggestions would be appreciated.

 For example (see attached script), first I create the DA ('dof' and
 'width' are defined):

 ?da = PETSc.DA().create(dim=2, sizes=[8, 8], dof=dof,
 ? ? ? ? ? ? ? ? ? ? ? ? boundary_type='periodic',
 ? ? ? ? ? ? ? ? ? ? ? ? stencil_width=width,
 ? ? ? ? ? ? ? ? ? ? ? ? stencil_type='box')

 Next, I create global and local vectors:

 ?gvec = da.createGlobalVec()
 ?lvec = da.createLocalVector()

 Then, put something into the local vectors, and scatter:

 ?with da.getVecArray(lvec) as va:
 ? ?if dof == 1:
 ? ? ?va.array[width:-width,width:-width] = MPI.COMM_WORLD.rank+1
 ? ?else:
 ? ? ?va.array[width:-width,width:-width,0] = MPI.COMM_WORLD.rank+1
 ?da.localToGlobal(lvec, gvec)

 Finally, get the local vector again and print:

 ?da.globalToLocal(gvec, lvec)
 ?with da.getVecArray(lvec) as va:
 ? ?if dof == 1:
 ? ? ?print va.array[:,:]
 ? ?else:
 ? ? ?print va.array[:,:,0]

 With 'dof=1' I get the following on the second processor:

 [[ 1. ?2. ?2. ?2. ?2. ?1.]
 ?[ 1. ?2. ?2. ?2. ?2. ?1.]
 ?[ 1. ?2. ?2. ?2. ?2. ?1.]
 ?[ 1. ?2. ?2. ?2. ?2. ?1.]
 ?[ 1. ?2. ?2. ?2. ?2. ?1.]
 ?[ 1. ?2. ?2. ?2. ?2. ?1.]
 ?[ 1. ?2. ?2. ?2. ?2. ?1.]
 ?[ 1. ?2. ?2. ?2. ?2. ?1.]
 ?[ 1. ?2. ?2. ?2. ?2. ?1.]
 ?[ 1. ?2. ?2. ?2. ?2. ?1.]]

 With 'dof=3' I get the following on the second processor:

 [[ 0. ?0. ?0. ?0. ?0. ?0.]
 ?[ 0. ?0. ?0. ?0. ?2. ?0.]
 ?[ 0. ?0. ?0. ?0. ?2. ?0.]
 ?[ 0. ?0. ?0. ?2. ?2. ?0.]
 ?[ 0. ?0. ?0. ?2. ?2. ?0.]
 ?[ 0. ?0. ?0. ?2. ?2. ?0.]
 ?[ 0. ?0. ?0. ?2. ?2. ?0.]
 ?[ 0. ?0. ?0. ?2. ?2. ?2.]
 ?[ 0. ?0. ?0. ?2. ?2. ?2.]
 ?[ 0. ?0. ?0. ?0. ?0. ?2.]]

 I was expecting to see the same as the 'dof=1' case, since I am
 setting (and printing) va.array[:,:,0].


 Again, any suggestions would be appreciated. ?Thanks,
 Matt

 PS, thanks to the developers of PETSc and petsc4py for all your hard
 work. ?I had the pleasure of meeting Jed and Matt at KAUST recently.


 Mmm, it seems that this functionality is broken for the case of
 boundary_type='periodic'. I'll take a look.


No, sorry... It has nothing to do with periodicity... It is actually a
C/Fortran ordering issue I need to fix...

In Fortran 90, it seems you index a DA Vec array as A[dof,x,y,z]... ,
However, I think that for Python we should follow a more C-ish
indexing A[x,y,z,dof]. Or we could do it like PETSc in C, that is
A[z,y,x,dof] (wich is the transpose of the Fortran way) but it is
counter-intuitive to C (and likely Python) programmers ...

What do others think about this?




-- 
Lisandro Dalcin
---
CIMEC (INTEC/CONICET-UNL)
Predio CONICET-Santa Fe
Colectora RN 168 Km 472, Paraje El Pozo
3000 Santa Fe, Argentina
Tel: +54-342-4511594 (ext 1011)
Tel/Fax: +54-342-4511169


[petsc-users] petsc4py DA with dof1 confusion

2011-07-04 Thread Lisandro Dalcin
On 4 July 2011 16:19, Matthew Emmett memmett at unc.edu wrote:
 On Mon, Jul 4, 2011 at 2:52 PM, Lisandro Dalcin dalcinl at gmail.com wrote:
 No, sorry... It has nothing to do with periodicity... It is actually a
 C/Fortran ordering issue I need to fix...

 Ah, I see. ?Thanks for looking into this so quickly. ?Let me know if
 there is anything I can do on my end to help. ?I will poke around the
 code for petsc4py, but you obviously know it better than I do.

 In Fortran 90, it seems you index a DA Vec array as A[dof,x,y,z]... ,
 However, I think that for Python we should follow a more C-ish
 indexing A[x,y,z,dof]. Or we could do it like PETSc in C, that is
 A[z,y,x,dof] (wich is the transpose of the Fortran way) but it is
 counter-intuitive to C (and likely Python) programmers ...

 What do others think about this?

 I think A[x,y,z,dof] is probably the most intuitive for Python.


OK. I pushed a fix solving the issue. As a consequence, the underlying
NumPy array do have a mixed ordering (no C nor Fortran) where the
latest dim varies fastest, then first, second, third. So, you have to
be very careful about accessing the underlying array (I mean, doing
va.array[...] in your previous code), I would recomend you against
that, and code like this instead:


  with da.getVecArray(lvec) as va:
if dof == 1:
  va[:,:] = rank+1
else:
  va[:,:,0] = rank+1
  da.localToGlobal(lvec, gvec)

  da.globalToLocal(gvec, lvec)
  with da.getVecArray(lvec) as va:
if dof == 1:
  print va[:,:]
else:
  print va[:,:,0]




-- 
Lisandro Dalcin
---
CIMEC (INTEC/CONICET-UNL)
Predio CONICET-Santa Fe
Colectora RN 168 Km 472, Paraje El Pozo
3000 Santa Fe, Argentina
Tel: +54-342-4511594 (ext 1011)
Tel/Fax: +54-342-4511169


[petsc-users] Flexiable AIJ matrix for nonzeros -- hash table version

2011-05-10 Thread Lisandro Dalcin
2011/5/9 Gong Ding gdiso at ustc.edu:
 Hi
 I created a flexible AIJ matrix type ?named as FAIJ. It has an extra hash 
 table for nonzero which do not have pre-allocated position. And the buffered 
 values in hash table will be flushed to AIJ array at MatAssemblyEnd. Both 
 sequential and parallel version have been implemented and tested.


Great!! Have you thought about extending regular AIJ matrices to use a
hash table on first assembly (perhaps activated with an option) and
switch back to regular assembly for the next assemblies? I'm assuming
here that calls to MatSetValues() fill the hash table even if values
are zero.

 The test results show that for middle size problem ? matrix size around 1M, 
 FAIJ without any pre-allocation still pretty fast. The shortcoming is memory 
 overkill exists. Fortunately, modern computers have larger and larger memory. 
 ?User with FAIJ matrix only needs to provide approximate nonzero pattern or 
 even skip this step.


How did you wrote your hash table? What are the keys and values? Are
you using a hash table per row? Or it is actuallyan  (i,j)-a hash
table?

 However, direct solvers such as mumps and superlu_dist will crash. This is 
 caused by type compare functions, i.e.:
 ierr = PetscTypeCompare((PetscObject)A,MATMPIAIJ,isAIJ);CHKERRQ(ierr);
 FAIJ matrix is derived from AIJ, but it has its own type name, i.e. 
 MATMPIFAIJ. The above function will assign isAIJ false.


That's the reason I think that adding this feature to regular AIJ
matrices could be better.

 Is there any function like ?PetscBaseType? exist? ?FAIJ matrix can hold its 
 type name as well as its base type name.

 I hope FAIJ can be accepted by petsc 3.2.
 Thanks.


I think Satish pushed something related to this a few days ago...


-- 
Lisandro Dalcin
---
CIMEC (INTEC/CONICET-UNL)
Predio CONICET-Santa Fe
Colectora RN 168 Km 472, Paraje El Pozo
3000 Santa Fe, Argentina
Tel: +54-342-4511594 (ext 1011)
Tel/Fax: +54-342-4511169


[petsc-users] Building with MKL 10.3

2011-03-14 Thread Lisandro Dalcin
On 14 March 2011 20:48, Barry Smith bsmith at mcs.anl.gov wrote:

 ?Eric,

 ? With the current PETSc design that uses one MPI process per core I am not 
 sure that it makes sense to use multi-threaded MKL since the processes are 
 already using all the cores. I could be wrong.

 ? Barry


Well, I think that PETSc should support sequential codes using
multi-threaded blas-lapack and perhaps some OpenMP. Moreover, perhaps
we could get some little speedup in MatMult for these kinds of
users/applications.

-- 
Lisandro Dalcin
---
CIMEC (INTEC/CONICET-UNL)
Predio CONICET-Santa Fe
Colectora RN 168 Km 472, Paraje El Pozo
3000 Santa Fe, Argentina
Tel: +54-342-4511594 (ext 1011)
Tel/Fax: +54-342-4511169


[petsc-users] Using PETSc from MATLAB code, experimental

2011-01-13 Thread Lisandro Dalcin
On 13 January 2011 08:38, Michel Cancelliere fernandez858 at gmail.com wrote:
 Hi Barry,
 is Matlab under Windows with petsc-cygwin supported?


Unlikely, we would need to build PETSc as a DLL... If anyone can
manage to do that, then Barry's work should work out of the box.

Or perhaps PETSc do build as a DLL under cygwin?

 Thanks,
 Michel
 On Sun, Dec 26, 2010 at 5:17 AM, Barry Smith bsmith at mcs.anl.gov wrote:

 ?PETSc users,

 ? It is now possible to write MATLAB programs (sequential) that use PETSc
 KSP, SNES, and TS solvers directly in MATLAB. The code is still experimental
 and incomplete. But if you are interested in trying it out, get the
 development release of PETSc
 http://www.mcs.anl.gov/petsc/petsc-as/developers/index.html join the
 development mailing list petsc-dev
 http://www.mcs.anl.gov/petsc/petsc-as/miscellaneous/mailing-lists.html, read
 bin/matlab/classes/PetscInitialize.m, configure and make PETSc and join the
 fun. We are definitely in need of more developers for this code.

 ? Barry







-- 
Lisandro Dalcin
---
CIMEC (INTEC/CONICET-UNL)
Predio CONICET-Santa Fe
Colectora RN 168 Km 472, Paraje El Pozo
Tel: +54-342-4511594 (ext 1011)
Tel/Fax: +54-342-4511169


[petsc-users] petsc rpath

2010-11-24 Thread Lisandro Dalcin
Perhaps you should add a --with-rpath option to ./configure ? For
reference, MPICH2 did that (=1.3), I believe the idea was to make the
life of distro packagers easier.

On 24 November 2010 15:28, Jens-Malte Gottfried jmgottfried at web.de wrote:
 Thanks a lot Satish, your patch did the job :-)

 Am 24.11.2010 17:19, schrieb Satish Balay:
 On Wed, 24 Nov 2010, Satish Balay wrote:

 B1;2601;0cOn Wed, 24 Nov 2010, Jed Brown wrote:

 On Wed, Nov 24, 2010 at 16:26, Jens-Malte Gottfried jmgottfried at 
 web.dewrote:

 I am trying to write a gentoo ebuild for petsc (my current work is
 already in the science overlay).
 Is there any possibility to disable writing the RPATH of the
 to-be-generated libpetsc.so?
 Currently, the library is built fine but RPATH is set to


 /usr/lib64:/usr/lib/gcc/x86_64-pc-linux-gnu/4.4.4:/usr/x86_64-pc-linux-gnu/lib
 which is not neccessary. Before trying to patch the build system,
 perhaps there is a simple configure option which I did not find.

 The standard build system does not provide an easy way to do this. ?If you
 are using petsc-dev, you use CMake (configure as usual, but then run make
 from $PETSC_DIR/$PETSC_ARCH instead of from $PETSC_DIR). ?The CMake build
 files default to setting RPATH, but it is easy to set the cache variables 
 to
 disable RPATH.
 Some of the above comes from 'gfortran -v' - which checkFortranLibraries() 
 etc do..

 For such known systems with known compatibility libraries [ and each
 language compiler can find each others libs] you can use:

 [for c++ build]
 I meant c build..

 -with-clib-autodetect=0 -with-fortranlib-autodetect=0 LIBS=-lgfortran

 [for c++ build]
 -with-clib-autodetect=0 -with-cxxlib-autodetect=0 
 -with-fortranlib-autodetect=0 LIBS='-lstdc++ -lgfortran'

 Wrt -rpath for -lpetsc, perhaps commenting out the following line in
 config/BuildSystem/config/setCompilers.py is sufficient:

 ? ? self.executeTest(self.checkSharedLinkerPaths)
 Actually you need '-L' there - so the following change will remove -rpath

 asterix:/home/balay/tmp/petsc-dist-test/config/BuildSystemhg diff
 diff -r 5fac76781491 config/setCompilers.py
 --- a/config/setCompilers.py ? ?Mon Sep 20 16:15:14 2010 -0500
 +++ b/config/setCompilers.py ? ?Wed Nov 24 10:12:46 2010 -0600
 @@ -1188,7 +1188,7 @@
 ? ? ? ?self.pushLanguage(language)
 ? ? ? ?# test '-R' before '-rpath' as sun compilers [c,fortran] don't give 
 proper errors with wrong options.
 ? ? ? ?if not Configure.isDarwin():
 - ? ? ? ?testFlags = ['-Wl,-rpath,', '-R','-rpath ' , '-Wl,-R,']
 + ? ? ? ?testFlags = []
 ? ? ? ?else:
 ? ? ? ? ?testFlags = []
 ? ? ? ?# test '-R' before '-Wl,-rpath' for SUN compilers [as cc on linux 
 accepts -Wl,-rpath, but ?f90  CC do not.
 asterix:/home/balay/tmp/petsc-dist-test/config/BuildSystem



 Satish





-- 
Lisandro Dalcin
---
CIMEC (INTEC/CONICET-UNL)
Predio CONICET-Santa Fe
Colectora RN 168 Km 472, Paraje El Pozo
Tel: +54-342-4511594 (ext 1011)
Tel/Fax: +54-342-4511169


[petsc-users] access da local vector data array in petsc4py

2010-10-12 Thread Lisandro Dalcin
On 12 October 2010 05:21, Amal Alghamdi amal.ghamdi at kaust.edu.sa wrote:
 Dear All,
 I have a question regarding petsc4py. I'd like to access the data array in a
 local vector generated from a one-dimension da, which has dof = 1. I used ?q
 = localVec.getArray(). However, q shape now does not take into account dof
 which is one in my example. I always need to reshape q to match the intended
 dof. Is there a way to do that automatically?
 Thank you very much.
 Amal

Almost all functions and methods in petsc4py treat arrays as linear
data, I mean, unidimensional arrays. As Vec.getArray() knows nothing
about the DA, there is no way to get the DA shape. Note that recently
I've added DA.getVecArray() to provide some support for global
indexing, this new thing is not finished, I want to improve
interoperability with regular numpy arrays, the functionality you want
could be added there.

Other possibility would be to add a method to Vec, let say
Vec.getArrayBlock(), that return an array reshaped to (n/bs, bs) where
n is the local vec size and bs is the block size (should equal
ndof for Vec's originating from DA). What do you think?

In the mean time, you could use an utility function like the one below:

def get_array_block(vec):
n = vec.getLocalSize()
bs = vec.getBlockSize()
a = numpy.asarray(vec)
a.shape = (n//bs, bs)
return a


-- 
Lisandro Dalcin
---
CIMEC (INTEC/CONICET-UNL)
Predio CONICET-Santa Fe
Colectora RN 168 Km 472, Paraje El Pozo
Tel: +54-342-4511594 (ext 1011)
Tel/Fax: +54-342-4511169


[petsc-users] petsc4py

2010-09-30 Thread Lisandro Dalcin
2010/9/29 ??? ??? kutuzovnp at gmail.com:
 Where can i get it? I've tried this one --
 http://code.google.com/p/petsc4py/source/browse/demo/ode/rober.py


This is from branch release-1.1, a  new release 1.1.2 is planned for
the weekend. As you can see from the output at the very end, SNES is
using matrix free. Note that I've not changed he code code, just
passed -snes_mf

What petsc4py and PETSc versions are you using?

[dalcinl at trantor petsc4py-release-1.1]$ python demo/ode/rober.py
-ts_view -snes_mf
TS Object:
  type: theta
Theta=0.5
Extrapolation=no
  maximum steps=100
  maximum time=1e+30
  total number of nonlinear solver iterations=203
  total number of linear solver iterations=311
  SNES Object:
type: ls
  line search variant: SNESLineSearchCubic
  alpha=0.0001, maxstep=1e+08, minlambda=1e-12
maximum iterations=50, maximum function evaluations=1
tolerances: relative=1e-08, absolute=1e-50, solution=1e-08
total number of linear solver iterations=3
total number of function evaluations=8
KSP Object:
  type: gmres
GMRES: restart=30, using Classical (unmodified) Gram-Schmidt
Orthogonalization with no iterative refinement
GMRES: happy breakdown tolerance 1e-30
  maximum iterations=1, initial guess is zero
  tolerances:  relative=1e-05, absolute=1e-50, divergence=1
  left preconditioning
  using PRECONDITIONED norm type for convergence test
PC Object:
  type: none
  linear system matrix = precond matrix:
  Matrix Object:
type=mffd, rows=3, cols=3
matrix-free approximation:
  err=1e-07 (relative error in function evaluation)
  Using wp compute h routine
  Computes normA
  Does not compute normU


-- 
Lisandro Dalcin
---
CIMEC (INTEC/CONICET-UNL)
Predio CONICET-Santa Fe
Colectora RN 168 Km 472, Paraje El Pozo
Tel: +54-342-4511594 (ext 1011)
Tel/Fax: +54-342-4511169


[petsc-users] petsc4py, plotting DA and writing to file

2010-09-27 Thread Lisandro Dalcin
On 27 September 2010 13:33, Barry Smith bsmith at mcs.anl.gov wrote:

 ?Lisandro

 ? ?In C he could do VecView(vec,PETSC_VIEWER_DRAW_WORLD); can he not do the 
 equivalent in Python?
 What syntax would he use?

viewer = PETSc.Viewer.DRAW(globalvec.comm)
globalvec.view(viewer) # or viewer(globalvec)


 Perhaps the 2d bratu python example could demonstrate this since it is simple 
 and there is no monkeying with copying to one process etc.


Well, that example does not even use a DA, the plot is going to be a bit ugly.

 ? Barry

 On Sep 27, 2010, at 11:27 AM, Lisandro Dalcin wrote:

 On 27 September 2010 13:12, Amal Alghamdi amal.ghamdi at kaust.edu.sa 
 wrote:
 Dear All,
 I'm using the DA structure in petsc4py. I'd like to know please what is the
 right way to:
 plot the DA vector (the global vector).
 write the global vector to a file.
 Is that each process writes or draws its own part? or I should communicate
 all the data to one process? or none of these?!!
 Thank you very much
 Amal

 Use a PETSc.Viewer().createBinary() to save the global vector. Each
 process will save its own part.

 In order to plot it, I think you should use numpy.fromfile() to load
 the data, and then plot the array.

 Other way would be to use a PETSc.Scatter.toZero() to get the data at
 process 0, and then plot it

 scatter, seqvec = PETSc.Scatter.toZero(globalvec)
 im = PETSc.InsertMode.INSERT_VALUES
 sm = PETSc.ScatterMode.FORWARD
 scatter.scatter(globalvec, seqvec, im, sm)
 if globalvec.comm.rank == 0:
 ? ?plot(seqvec.array)


 --
 Lisandro Dalcin
 ---
 CIMEC (INTEC/CONICET-UNL)
 Predio CONICET-Santa Fe
 Colectora RN 168 Km 472, Paraje El Pozo
 Tel: +54-342-4511594 (ext 1011)
 Tel/Fax: +54-342-4511169





-- 
Lisandro Dalcin
---
CIMEC (INTEC/CONICET-UNL)
Predio CONICET-Santa Fe
Colectora RN 168 Km 472, Paraje El Pozo
Tel: +54-342-4511594 (ext 1011)
Tel/Fax: +54-342-4511169


[petsc-users] petsc4py, plotting DA and writing to file

2010-09-27 Thread Lisandro Dalcin
On 27 September 2010 14:02, Barry Smith bsmith at mcs.anl.gov wrote:

 On Sep 27, 2010, at 11:47 AM, Lisandro Dalcin wrote:

 On 27 September 2010 13:33, Barry Smith bsmith at mcs.anl.gov wrote:

 ?Lisandro

 ? ?In C he could do VecView(vec,PETSC_VIEWER_DRAW_WORLD); can he not do the 
 equivalent in Python?
 What syntax would he use?

 viewer = PETSc.Viewer.DRAW(globalvec.comm)
 globalvec.view(viewer) # or viewer(globalvec)


 Perhaps the 2d bratu python example could demonstrate this since it is 
 simple and there is no monkeying with copying to one process etc.


 Well, that example does not even use a DA, the plot is going to be a bit 
 ugly.

 ?Sorry. Barry


I could add a poisson2d implemented with DA and working in parallel.



 ? Barry

 On Sep 27, 2010, at 11:27 AM, Lisandro Dalcin wrote:

 On 27 September 2010 13:12, Amal Alghamdi amal.ghamdi at kaust.edu.sa 
 wrote:
 Dear All,
 I'm using the DA structure in petsc4py. I'd like to know please what is 
 the
 right way to:
 plot the DA vector (the global vector).
 write the global vector to a file.
 Is that each process writes or draws its own part? or I should communicate
 all the data to one process? or none of these?!!
 Thank you very much
 Amal

 Use a PETSc.Viewer().createBinary() to save the global vector. Each
 process will save its own part.

 In order to plot it, I think you should use numpy.fromfile() to load
 the data, and then plot the array.

 Other way would be to use a PETSc.Scatter.toZero() to get the data at
 process 0, and then plot it

 scatter, seqvec = PETSc.Scatter.toZero(globalvec)
 im = PETSc.InsertMode.INSERT_VALUES
 sm = PETSc.ScatterMode.FORWARD
 scatter.scatter(globalvec, seqvec, im, sm)
 if globalvec.comm.rank == 0:
 ? ?plot(seqvec.array)


 --
 Lisandro Dalcin
 ---
 CIMEC (INTEC/CONICET-UNL)
 Predio CONICET-Santa Fe
 Colectora RN 168 Km 472, Paraje El Pozo
 Tel: +54-342-4511594 (ext 1011)
 Tel/Fax: +54-342-4511169





 --
 Lisandro Dalcin
 ---
 CIMEC (INTEC/CONICET-UNL)
 Predio CONICET-Santa Fe
 Colectora RN 168 Km 472, Paraje El Pozo
 Tel: +54-342-4511594 (ext 1011)
 Tel/Fax: +54-342-4511169





-- 
Lisandro Dalcin
---
CIMEC (INTEC/CONICET-UNL)
Predio CONICET-Santa Fe
Colectora RN 168 Km 472, Paraje El Pozo
Tel: +54-342-4511594 (ext 1011)
Tel/Fax: +54-342-4511169


[petsc-users] petsc4py, plotting DA and writing to file

2010-09-27 Thread Lisandro Dalcin
On 27 September 2010 13:54, Jed Brown jed at 59a2.org wrote:
 On Mon, Sep 27, 2010 at 18:47, Lisandro Dalcin dalcinl at gmail.com wrote:
 viewer = PETSc.Viewer.DRAW(globalvec.comm)
 globalvec.view(viewer) # or viewer(globalvec)

 What do you think about having a top-level matplotlib viewer?

Yes, we could have it, and a VisIt plugin, and a ParaView one, and
what about Mayavi, and put here your preferred viewer... This could
go to a separate petsc4py.plotting module (I prefer to not contaminate
the petsc4py.PETSc namespace with stuff that is not pure-PETSc)

?Ideally
 it would have an option to drop you into an interactive python session
 after the initial view.


Could you elaborate?

 I know this brings up the multiple dispatch problem, since the viewer
 would naturally be distributed with petsc4p, but then you don't get to
 modify the case statement in VecView_XX.


Sorry, now I'm confused. Are you suggesting this for petsc4py, or core PETSc?

-- 
Lisandro Dalcin
---
CIMEC (INTEC/CONICET-UNL)
Predio CONICET-Santa Fe
Colectora RN 168 Km 472, Paraje El Pozo
Tel: +54-342-4511594 (ext 1011)
Tel/Fax: +54-342-4511169


[petsc-users] petsc4py, plotting DA and writing to file

2010-09-27 Thread Lisandro Dalcin
On 27 September 2010 14:34, Amal Alghamdi amal.ghamdi at kaust.edu.sa wrote:
 Thank you very much everybody for the prompt help.
 Amal


I've just pushed a new demo: matrix-free Poisson 2D

http://code.google.com/p/petsc4py/source/browse/demo/poisson2d/poisson2d.py


At the very end you have:

draw = PETSc.Viewer.DRAW(x.comm)
OptDB['draw_pause'] = 1
draw(x)

Additionally, I've updated other demos. Take a look here for
matplotlib usage gathering values to processor 0:

http://code.google.com/p/petsc4py/source/browse/demo/bratu3d/bratu3d.py

(note that this demo is 3D, I'm plotting the solution at plane Z=0.5)



-- 
Lisandro Dalcin
---
CIMEC (INTEC/CONICET-UNL)
Predio CONICET-Santa Fe
Colectora RN 168 Km 472, Paraje El Pozo
Tel: +54-342-4511594 (ext 1011)
Tel/Fax: +54-342-4511169


  1   2   >