Re: [petsc-users] Periodic BC in petsc DMPlex / firedrake

[Matthew Knepley]

On Wed, Oct 24, 2018 at 11:36 AM Maximilian Hartig 
wrote:

>
>
> On 24. Oct 2018, at 12:49, Matthew Knepley  wrote:
>
> On Wed, Oct 24, 2018 at 6:29 AM Lawrence Mitchell  wrote:
>
>> Hi Max,
>>
>> (I'm cc'ing in the petsc-users mailing list which may have more advice,
>> if you are using PETSc you should definitely subscribe!
>>
>> > On 24 Oct 2018, at 09:27, Maximilian Hartig 
>> wrote:
>> >
>> > Hello Lawrence,
>> >
>> > sorry to message you out of the blue. My name is Max and I found your
>> post on GitHub (https://github.com/firedrakeproject/firedrake/issues/1246
>>  ) on DMPlex being able to read periodic gmsh files. I am currently
>> trying to do just that (creating a periodic DMPlex mesh with gmsh) in the
>> context of my PhD work. So far I haven’ t found any documentation on the
>> periodic BC’s with DMPlex and gmsh in the official petsc documentation.
>> > I was wondering whether you’d be so kind as to point me in a general
>> direction concerning how to achieve this. You seem experienced in using
>> petsc and I would greatly appreciate your help.
>>
>>
>> I think the answer is "it depends". If you're just using DMPlex directly
>> and all the of the functionality with PetscDS, then I /think/ that reading
>> periodic meshes via gmsh (assuming you're using the appropriate gmsh mesh
>> format [v2]) "just works".
>>
>
> There are two phases here: topological and geometric. DMPlex represents
> the periodic topological entity directly. For example,  a circle is just a
> segment with one end hooked to the other. Vertices are not duplicated, or
> mapped to each other. This makes topology simple and easy to implement.
> However, then geometry is more complicated. What Plex does is allow
> coordinates to be represented by a discontinuous field taking values on
> cells, in addition to vertices. In our circle example, each cells near the
> cut will have 2 coordinates, one for each vertex, but they will not agree
> across the cut. If you define a periodic domain, then Plex can construct
> this coordinate field automatically using DMPlexLocalize(). These DG
> coordinates are then used by the integration routines.
>
>
> Ok, I think I understand the concept. DMPlex reads information about both
> topology and coordinates from the .msh file. Creating a periodic mesh in
> gmsh then should allow DMPlex to identify the periodic boundaries as the
> “cut” and build the mesh topology accordingly. Coordinate information is
> handled separately.
> That means, as Lawrence suggested, building periodic meshes in gmsh and
> reading them in to petsc’s DMPlex should indeed “just work”.  (From the
> user perspective). The only extra step is to call DMLocalizeCoordinates()
> after DMPlexReadFromFile(). Sorry to reiterate, I am just trying to make
> sense of this.
>
>
>
>> From my side, the issue is to do with mapping that coordinate field into
>> one that I understand (within Firedrake). You may not have this problem.
>>
>
> Firedrake uses its own coordinate mapping and integration routines, so
> they must manage the second part independently. I hope to get change this
> slightly soon by making the Firedrake representation a DMField, so that it
> looks the same to Plex.
>
>   Thanks,
>
> Matt
>
>
>> Thanks,
>>
>> Lawrence
>
>
>
> --
> 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
>
> https://www.cse.buffalo.edu/~knepley/
>
>
> To read periodic meshes from GMSH, you need to use the option
> -dm_plex_gmsh_periodic and DMPlexCreateFromFile
>
>
> Ahh, thanks. I was missing the option " -dm_plex_gmsh_periodic “. But
> using this option I now generate a segmentation fault error when calling
> VecView() on the solution vector with vtk and hdf5 viewers. Any suggestions?
>

 Small example? VTK is deprecated. HDF5 should work, although it will
require you to have proper coordinates I think. We have to
think about what you mean. If its for a checkpoint, there is no problem,
but for viz, those programs do not understand periodicity. Thus I embed it
in a higher dimensional space.

   Matt

> See  src/dm/impls/plex/examples/tests/ex1.c. An example runs
>
> $ ./ex1 -filename
> ${PETSC_DIR}/share/petsc/datafiles/meshes/cube_periodic_bin.msh
> -dm_plex_gmsh_periodic -dm_view ::ascii_info_detail -interpolate -test_shape
>
> generating periodic meshes in gmsh may be tricky, Lisandro for sure may
> advice.
>
>
> Thanks,
> Max
>
>

-- 
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

https://www.cse.buffalo.edu/~knepley/ 

Re: [petsc-users] Using BDDC preconditioner for assembled matrices

[Abdullah Ali Sivas]

Hi Stefano,

I am trying to solve the div-div problem (or grad-div problem in strong
form) with a H(div)-conforming FEM. I am getting the matrices from an
external source (to be clear, from an ngsolve script) and I am not sure if
it is possible to get a MATIS matrix out of that. So I am just treating it
as if I am not able to access the assembly code. The results are 2, 31, 26,
27, 31 iterations, respectively, for matrix sizes 282, 1095, 4314, 17133,
67242, 267549. However, norm of the residual also grows significantly;
7.38369e-09 for 1095 and 5.63828e-07 for 267549. I can try larger sizes, or
maybe this is expected for this case.

As a side question, if we are dividing the domain into number of MPI
processes subdomains, does it mean that convergence is affected negatively
by the increasing number of processes? I know that alternating Schwarz
method and some other domain decomposition methods sometimes suffer from
the decreasing radius of the subdomains. It sounds like BDDC is pretty
similar to those by your description.

Best wishes,
Abdullah Ali Sivas

On Wed, 24 Oct 2018 at 05:28, Stefano Zampini 
wrote:

> Abdullah,
>
> The "Neumann" problems Jed is referring to result from assembling your
> problem on each subdomain ( = MPI process) separately.
> Assuming you are using FEM, these problems have been historically  named
> "Neumann" as they correspond to a problem with natural boundary conditions
> (Neumann bc for Poisson).
> Note that in PETSc the subdomain decomposition is associated with the mesh
> decomposition.
>
> When converting from an assembled AIJ matrix to a MATIS format, such
> "Neumann" information is lost.
> You can disassemble an AIJ matrix, in the sense that you can find local
> matrices A_j such that A = \sum_j R^T_j A_j R_j (as it is done in ex72.c),
> but you cannot guarantee (unless if you solve an optimization problem) that
> the disassembling will produce subdomain Neumann problems that are
> consistent with your FEM problem.
>
> I have added such disassembling code a few months ago, just to have
> another alternative for preconditioning AIJ matrices in PETSc; there are
> few tweaks one can do to improve the quality of the disassembling, but I
> discourage its usage unless you don't have access to the FEM assembly code.
>
> With that said, what problem are you trying to solve? Are you using DMDA
> or DMPlex? What are the results you obtained with using the automatic
> disassembling?
>
> Il giorno mer 24 ott 2018 alle ore 08:14 Abdullah Ali Sivas <
> abdullahasi...@gmail.com> ha scritto:
>
>> Hi Jed,
>>
>> Thanks for your reply. The assembled matrix I have corresponds to the
>> full problem on the full mesh. There are no "Neumann" problems (or any sort
>> of domain decomposition) defined in the code generates the matrix. However,
>> I think assembling the full problem is equivalent to implicitly assembling
>> the "Neumann" problems, since the system can be partitioned as;
>>
>> [A_{LL} | A_{LI}]  [u_L] [F]
>> ---|  = -
>> [A_{IL}  |A_{II} ]   [u_I]  [G]
>>
>> and G should correspond to the Neumann problem. I might be thinking wrong
>> (or maybe I completely misunderstood the idea), if so please correct me.
>> But I think that the problem is that I am not explicitly telling PCBDDC
>> which dofs are interface dofs.
>>
>> Regards,
>> Abdullah Ali Sivas
>>
>> On Tue, 23 Oct 2018 at 23:16, Jed Brown  wrote:
>>
>>> Did you assemble "Neumann" problems that are compatible with your
>>> definition of interior/interface degrees of freedom?
>>>
>>> Abdullah Ali Sivas  writes:
>>>
>>> > Dear all,
>>> >
>>> > I have a series of linear systems coming from a PDE for which BDDC is
>>> an
>>> > optimal preconditioner. These linear systems are assembled and I read
>>> them
>>> > from a file, then convert into MATIS as required (as in
>>> >
>>> https://www.mcs.anl.gov/petsc/petsc-current/src/ksp/ksp/examples/tutorials/ex72.c.html
>>> > ). I expect each of the systems converge to the solution in almost same
>>> > number of iterations but I don't observe it. I think it is because I
>>> do not
>>> > provide enough information to the preconditioner. I can get a list of
>>> inner
>>> > dofs and interface dofs. However, I do not know how to use them. Has
>>> anyone
>>> > have any insights about it or done something similar?
>>> >
>>> > Best wishes,
>>> > Abdullah Ali Sivas
>>>
>>
>
> --
> Stefano
>

Re: [petsc-users] Usage of MATMPISBAIJ

[Zhang, Hong]

Ian:
I am currently thinking of doing a Sparse-Dense matrix matrix multiplication, 
where the Sparse matrix is a larger symmetric matrix, and the dense matrix is a 
skinny matrix. The sparse matrix is used mainly because of the memory issue, 
and the sparsity is just around 10~20%. Thus I am thinking of using MPISBAIJ 
type to reduce the memory cost. Does MPISBAIJ support MatMatMult with MPIDense 
matrix? If so, how would the data be distributed over processors? Would it be 
the same as MPIAIJ (row are by default evenly distributed over processors, but 
that does not make sense to MPISBAIJ as I suppose the lower part of the matrix 
is not stored)?
We do not support MatMatMult with MPISBAIJ and MPIDense matrix, only MPIAIJ and 
MPIDense matrix.

Also, just another quick question, in your experience, considering the overhead 
cost for sparse algorithm, how sparse should a matrix be for the sparse 
algorithm to outperform dense algorithm in terms of computation time?
Sparse matrix normally is defined as the sparsity less than 10%, or less than 
1%. I would consider 20% as a dense matrix.
The complexity of C=A*B is O(nnzA*nnzB), nnzA = num of nonzeros of A.
Hong

Re: [petsc-users] PetscInt overflow

[Matthew Knepley]

As it says, if you are looking at performance, you should configure using
--with-debugging=0.

It says SLEPc is using 8GB. Is this what you see?

  Matt

On Wed, Oct 24, 2018 at 12:30 PM Jan Grießer 
wrote:

> I also run it with the -log_summary :
> -- PETSc Performance Summary:
> --
>
>
>
>   ##
>   ##
>   #  WARNING!!!#
>   ##
>   #   This code was compiled with a debugging option,  #
>   #   To get timing results run ./configure#
>   #   using --with-debugging=no, the performance will  #
>   #   be generally two or three times faster.  #
>   ##
>   ##
>
>
> /work/ws/nemo/fr_jg1080-FreeSurface_Glass-0/GlassSystems/PeriodicSystems/N50T0.001/SolveEigenvalueProblem_par/Test/Eigensolver_petsc_slepc_no_argparse.py
> on a arch-linux2-c-debug named int02.nemo.privat with 20 processors, by
> fr_jg1080 Wed Oct 24 18:26:30 2018
> Using Petsc Release Version 3.9.4, Sep, 11, 2018
>
>  Max   Max/MinAvg  Total
> Time (sec):   7.474e+02  1.0   7.474e+02
> Objects:  3.600e+01  1.0   3.600e+01
> Flop: 1.090e+11  1.00346   1.089e+11  2.177e+12
> Flop/sec:1.459e+08  1.00346   1.457e+08  2.913e+09
> Memory:   3.950e+08  1.00296  7.891e+09
> MPI Messages: 3.808e+04  1.0   3.808e+04  7.615e+05
> MPI Message Lengths:  2.211e+10  1.00217   5.802e+05  4.419e+11
> MPI Reductions:   1.023e+05  1.0
>
> Flop counting convention: 1 flop = 1 real number operation of type
> (multiply/divide/add/subtract)
> e.g., VecAXPY() for real vectors of length N
> --> 2N flop
> and VecAXPY() for complex vectors of length N
> --> 8N flop
>
> Summary of Stages:   - Time --  - Flop -  --- Messages
> ---  -- Message Lengths --  -- Reductions --
> Avg %Total Avg %Total   counts
>  %Total Avg %Total   counts   %Total
>  0:  Main Stage: 7.4739e+02 100.0%  2.1773e+12 100.0%  7.615e+05
> 100.0%  5.802e+05  100.0%  1.022e+05 100.0%
>
>
> 
> See the 'Profiling' chapter of the users' manual for details on
> interpreting output.
> Phase summary info:
>Count: number of times phase was executed
>Time and Flop: Max - maximum over all processors
>Ratio - ratio of maximum to minimum over all processors
>Mess: number of messages sent
>Avg. len: average message length (bytes)
>Reduct: number of global reductions
>Global: entire computation
>Stage: stages of a computation. Set stages with PetscLogStagePush() and
> PetscLogStagePop().
>   %T - percent time in this phase %F - percent flop in this
> phase
>   %M - percent messages in this phase %L - percent message lengths
> in this phase
>   %R - percent reductions in this phase
>Total Mflop/s: 10e-6 * (sum of flop over all processors)/(max time over
> all processors)
>
> 
>
>
>   ##
>   ##
>   #  WARNING!!!#
>   ##
>   #   This code was compiled with a debugging option,  #
>   #   To get timing results run ./configure#
>   #   using --with-debugging=no, the performance will  #
>   #   be generally two or three times faster.  #
>   ##
>   ##
>
>
> EventCount  Time (sec) Flop
>  --- Global ---  --- Stage ---   Total
>Max Ratio  Max Ratio   Max  Ratio  Mess   Avg len
> Reduct  %T %F %M %L %R  %T %F %M %L %R Mflop/s
>
> 
>
> --- Event Stage 0: Main Stage
>
> BuildTwoSidedF 2 1.0 2.6670e-01 1.0 0.00e+00 0.0 0.0e+00 0.0e+00
> 0.0e+00  0  0  0  0  0   0  0  0  0  0 0
> VecSet 2 1.0 6.8650e-03 1.8 0.00e+00 0.0 0.0e+00 0.0e+00
> 0.0e+00  0  0  0  0  0  

Re: [petsc-users] PetscInt overflow

[Jan Grießer]

I also run it with the -log_summary :
-- PETSc Performance Summary:
--



  ##
  ##
  #  WARNING!!!#
  ##
  #   This code was compiled with a debugging option,  #
  #   To get timing results run ./configure#
  #   using --with-debugging=no, the performance will  #
  #   be generally two or three times faster.  #
  ##
  ##


/work/ws/nemo/fr_jg1080-FreeSurface_Glass-0/GlassSystems/PeriodicSystems/N50T0.001/SolveEigenvalueProblem_par/Test/Eigensolver_petsc_slepc_no_argparse.py
on a arch-linux2-c-debug named int02.nemo.privat with 20 processors, by
fr_jg1080 Wed Oct 24 18:26:30 2018
Using Petsc Release Version 3.9.4, Sep, 11, 2018

 Max   Max/MinAvg  Total
Time (sec):   7.474e+02  1.0   7.474e+02
Objects:  3.600e+01  1.0   3.600e+01
Flop: 1.090e+11  1.00346   1.089e+11  2.177e+12
Flop/sec:1.459e+08  1.00346   1.457e+08  2.913e+09
Memory:   3.950e+08  1.00296  7.891e+09
MPI Messages: 3.808e+04  1.0   3.808e+04  7.615e+05
MPI Message Lengths:  2.211e+10  1.00217   5.802e+05  4.419e+11
MPI Reductions:   1.023e+05  1.0

Flop counting convention: 1 flop = 1 real number operation of type
(multiply/divide/add/subtract)
e.g., VecAXPY() for real vectors of length N
--> 2N flop
and VecAXPY() for complex vectors of length N
--> 8N flop

Summary of Stages:   - Time --  - Flop -  --- Messages ---
-- Message Lengths --  -- Reductions --
Avg %Total Avg %Total   counts
 %Total Avg %Total   counts   %Total
 0:  Main Stage: 7.4739e+02 100.0%  2.1773e+12 100.0%  7.615e+05
100.0%  5.802e+05  100.0%  1.022e+05 100.0%


See the 'Profiling' chapter of the users' manual for details on
interpreting output.
Phase summary info:
   Count: number of times phase was executed
   Time and Flop: Max - maximum over all processors
   Ratio - ratio of maximum to minimum over all processors
   Mess: number of messages sent
   Avg. len: average message length (bytes)
   Reduct: number of global reductions
   Global: entire computation
   Stage: stages of a computation. Set stages with PetscLogStagePush() and
PetscLogStagePop().
  %T - percent time in this phase %F - percent flop in this
phase
  %M - percent messages in this phase %L - percent message lengths
in this phase
  %R - percent reductions in this phase
   Total Mflop/s: 10e-6 * (sum of flop over all processors)/(max time over
all processors)



  ##
  ##
  #  WARNING!!!#
  ##
  #   This code was compiled with a debugging option,  #
  #   To get timing results run ./configure#
  #   using --with-debugging=no, the performance will  #
  #   be generally two or three times faster.  #
  ##
  ##


EventCount  Time (sec) Flop
 --- Global ---  --- Stage ---   Total
   Max Ratio  Max Ratio   Max  Ratio  Mess   Avg len
Reduct  %T %F %M %L %R  %T %F %M %L %R Mflop/s


--- Event Stage 0: Main Stage

BuildTwoSidedF 2 1.0 2.6670e-01 1.0 0.00e+00 0.0 0.0e+00 0.0e+00
0.0e+00  0  0  0  0  0   0  0  0  0  0 0
VecSet 2 1.0 6.8650e-03 1.8 0.00e+00 0.0 0.0e+00 0.0e+00
0.0e+00  0  0  0  0  0   0  0  0  0  0 0
VecScatterBegin 2002 1.0 1.4380e+01 1.0 0.00e+00 0.0 7.6e+05 5.8e+05
0.0e+00  2  0100100  0   2  0100100  0 0
VecScatterEnd   2002 1.0 3.7604e+01 1.5 0.00e+00 0.0 0.0e+00 0.0e+00
0.0e+00  4  0  0  0  0   4  0  0  0  0 0
VecSetRandom   1 1.0 1.6440e-01 1.0 0.00e+00 0.0 0.0e+00 0.0e+00
0.0e+00  0  0  0  0  0   0  0  0  0  0 0
MatMult  

[petsc-users] Usage of MATMPISBAIJ

[Ian Lin]

Hi PETSc team,

I am currently thinking of doing a Sparse-Dense matrix matrix multiplication, 
where the Sparse matrix is a larger symmetric matrix, and the dense matrix is a 
skinny matrix. The sparse matrix is used mainly because of the memory issue, 
and the sparsity is just around 10~20%. Thus I am thinking of using MPISBAIJ 
type to reduce the memory cost. Does MPISBAIJ support MatMatMult with MPIDense 
matrix? If so, how would the data be distributed over processors? Would it be 
the same as MPIAIJ (row are by default evenly distributed over processors, but 
that does not make sense to MPISBAIJ as I suppose the lower part of the matrix 
is not stored)? Also, just another quick question, in your experience, 
considering the overhead cost for sparse algorithm, how sparse should a matrix 
be for the sparse algorithm to outperform dense algorithm in terms of 
computation time?

Thanks a lot for your help.

Best regards,
Ian

Re: [petsc-users] PetscInt overflow

[Jan Grießer]

For some reason i get only this error message, also when is use the
-eps_view_pre. I started the program with nev=1, there the output is  (with
-bv_type vecs -bv_type mat -eps_view_pre)
EPS Object: 20 MPI processes
  type: krylovschur
50% of basis vectors kept after restart
using the locking variant
  problem type: symmetric eigenvalue problem
  selected portion of the spectrum: smallest real parts
  number of eigenvalues (nev): 1
  number of column vectors (ncv): 3
  maximum dimension of projected problem (mpd): 2
  maximum number of iterations: 1000
  tolerance: 1e-08
  convergence test: relative to the eigenvalue
BV Object: 20 MPI processes
  type: mat
  4 columns of global length 150
  vector orthogonalization method: classical Gram-Schmidt
  orthogonalization refinement: if needed (eta: 0.7071)
  block orthogonalization method: GS
  doing matmult as a single matrix-matrix product
DS Object: 20 MPI processes
  type: hep
  parallel operation mode: REDUNDANT
  solving the problem with: Implicit QR method (_steqr)
ST Object: 20 MPI processes
  type: shift
  shift: 0.
  number of matrices: 1




Am Mi., 24. Okt. 2018 um 16:14 Uhr schrieb Matthew Knepley <
knep...@gmail.com>:

> On Wed, Oct 24, 2018 at 10:03 AM Jan Grießer 
> wrote:
>
>> This is the error message i get from my program:
>> Shape of the dynamical matrix:  (150, 150)
>>
>
> Petsc installs a signal handler, so there should be a PETSc-specific
> message before this one. Is something eating
> your output?
>
>   Thanks,
>
>  Matt
>
>
>>
>> ===
>> =   BAD TERMINATION OF ONE OF YOUR APPLICATION PROCESSES
>> =   PID 122676 RUNNING AT n3512.nemo.privat
>> =   EXIT CODE: 9
>> =   CLEANING UP REMAINING PROCESSES
>> =   YOU CAN IGNORE THE BELOW CLEANUP MESSAGES
>>
>> ===
>>Intel(R) MPI Library troubleshooting guide:
>>   https://software.intel.com/node/561764
>>
>> ===
>>
>>
>> Am Mi., 24. Okt. 2018 um 16:01 Uhr schrieb Matthew Knepley <
>> knep...@gmail.com>:
>>
>>> On Wed, Oct 24, 2018 at 9:38 AM Jan Grießer 
>>> wrote:
>>>
 Hey,
 i tried to run my program as you said with -bv_type vecs and/or
 -bv_type mat, but instead of the PETScInt overflow i now get an MPI Error 9

>>>
>>> Send the actual error.
>>>
>>>   Thanks,
>>>
>>> Matt
>>>
>>>
 message, which i assume (after googling a little bit around) should be
 a memory problem. I tried to run it also on slightly bigger compute nodes
 on our cluster with 20 cores with each 12 GB and 24 GB but the problem
 still remains. The actual limit appears to be a 1.5 Million x 1.5 Million
 where i  searched for NEV = 1500 and MPD = 500/ 200 eigenvalues.
 Do you have maybe an idea what the error could be? I appended also the
 python method i used to solve the problem. I also tried to solve the
 problem with spectrum solving but the error message remains the same.

 def solve_eigensystem(DynMatrix_nn, NEV, MPD, Dimension):
 # Create the solver
 # E is used as an acronym for the EPS solver (EPS = Eigenvalue problem
 solver)
 E = SLEPc.EPS().create()

 # Set the preconditioner
 pc = PETSc.PC().create()
 pc.setType(pc.Type.BJACOBI)

 # Set the linear solver
 # Create the KSP object
 ksp = PETSc.KSP().create()
 # Create the solver, in this case GMRES
 ksp.setType(ksp.Type.GMRES)
 # Set the tolerances of the GMRES solver
 # Link it to the PC
 ksp.setPC(pc)

 # Set up the spectral transformations
 st = SLEPc.ST().create()
 st.setType("shift")
 st.setKSP(ksp)
 # MPD stands for "maximum projected dimension". It has to due with
 computational cost, please read Chap. 2.6.5 of SLEPc docu for
 # an explanation. At the moment mpd is only a guess
 E.setDimensions(nev=NEV, mpd = MPD)
 # Eigenvalues should be real, therefore we start to order them from the
 smallest real value |l.real|
 E.setWhichEigenpairs(E.Which.SMALLEST_REAL)
 # Since the dynamical matrix is symmetric and real it is hermitian
 E.setProblemType(SLEPc.EPS.ProblemType.HEP)
 # Use the Krylov Schur method to solve the eigenvalue problem
 E.setType(E.Type.KRYLOVSCHUR)
 # Set the convergence criterion to relative to the eigenvalue and the
 maximal number of iterations
 E.setConvergenceTest(E.Conv.REL)
 E.setTolerances(tol = 1e-8, max_it = 5000)
 # Set the matrix in order to solve
 E.setOperators(DynMatrix_nn, None)
 # Sets EPS options from the options database. This routine must be
 called before `setUp()` if the user is to be allowed to set dthe solver
 type.
 E.setFromOptions()
 # Sets up all the internal data structures necessary for the execution
 of 

Re: [petsc-users] Problems about Compiling Multifile Program

[Yingjie Wu]

Sorry, I made a mistake. I have a.C file with the same name in the current
directory. When I remove it, the compiler can be passed.
It seems that ${CLINKER} can automatically identify the extension of files
and compile them with different compilers.

Thanks for helping me.
Yingjie

Balay, Satish  于2018年10月24日周三 下午11:35写道：

> Please send complete logs. [i.e the complete compile output from your
> makefile. Also send your makefile]
>
> If your c++ code uses .cxx - it should get compiled with a c++ compiler.
>
> You might benefit by change CLINKER to CXXLINKER
>
> Satish
>
> On Wed, 24 Oct 2018, Yingjie Wu wrote:
>
> > Thank you very much for your reply.
> >
> > Because the files I added are written in c++, including definitions of
> > classes and functions. I encountered a mistake when I used your method.
> > myprogram: class1.o myprogram.o
> >${CLINKER} -o myprogram class1.o myprogram.o ${PETSC_LIB}
> > Error information is : unknown type name 'class'
> > It seems that I did not compile using the c++ compiler.
> >
> > How should I set it up so that it can compile successfully?
> >
> > Thanks,
> > Yingjie
> >
> > Matthew Knepley  于2018年10月24日周三 上午1:46写道：
> >
> > > On Tue, Oct 23, 2018 at 11:37 AM Yingjie Wu  wrote:
> > >
> > >> Dear Petsc developer:
> > >> Hi,
> > >> Thank you very much for your continuous help, I recently encountered
> some
> > >> difficulties in developing programs on Petsc.
> > >>
> > >> 1. I want to use my class definition (class1. h) and class functions
> > >> (class1. cpp) files in my Petsc program (myprogram. c) and compile my
> > >> program. I have written other program before:
> > >> g++ -c class1.cpp myprogram.c
> > >> g++ -o myprogram *.o
> > >> I also have some knowledge of Makefile, but I don't know how to
> modify it
> > >> to achieve compilation. I really want your help.
> > >>
> > >
> > > With PETSc makefiles, we have automatic rules, so you can use
> > >
> > > myprogram: class1.o myprogram.o
> > >${CLINKER} -o myprogram class1.o myprogram.o ${PETSC_LIB}
> > >
> > >
> > >> 2.  In my class definition and function implementation program, I
> want to
> > >> use Petsc data structures, such as PetscInt, PetscReal. In order to
> > >> maintain consistency of the main program. What kind of header files
> should
> > >> I add to my code, or do I need to add any code?
> > >>
> > >
> > > The easiest thing to do is just
> > >
> > > #include 
> > >
> > >   Thanks,
> > >
> > > Matt
> > >
> > >
> > >> Thanks,
> > >> Yingjie
> > >>
> > >
> > >
> > > --
> > > 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
> > >
> > > https://www.cse.buffalo.edu/~knepley/
> > > 
> > >
> >
>

Re: [petsc-users] Periodic BC in petsc DMPlex / firedrake

[Maximilian Hartig]

> On 24. Oct 2018, at 12:49, Matthew Knepley  wrote:
> 
> On Wed, Oct 24, 2018 at 6:29 AM Lawrence Mitchell  > wrote:
> Hi Max,
> 
> (I'm cc'ing in the petsc-users mailing list which may have more advice, if 
> you are using PETSc you should definitely subscribe!
> 
> > On 24 Oct 2018, at 09:27, Maximilian Hartig  > > wrote:
> > 
> > Hello Lawrence,
> > 
> > sorry to message you out of the blue. My name is Max and I found your post 
> > on GitHub (https://github.com/firedrakeproject/firedrake/issues/1246 
> >  ) on DMPlex 
> > being able to read periodic gmsh files. I am currently trying to do just 
> > that (creating a periodic DMPlex mesh with gmsh) in the context of my PhD 
> > work. So far I haven’ t found any documentation on the periodic BC’s with 
> > DMPlex and gmsh in the official petsc documentation. 
> > I was wondering whether you’d be so kind as to point me in a general 
> > direction concerning how to achieve this. You seem experienced in using 
> > petsc and I would greatly appreciate your help. 
> 
> 
> I think the answer is "it depends". If you're just using DMPlex directly and 
> all the of the functionality with PetscDS, then I /think/ that reading 
> periodic meshes via gmsh (assuming you're using the appropriate gmsh mesh 
> format [v2]) "just works".
> 
> There are two phases here: topological and geometric. DMPlex represents the 
> periodic topological entity directly. For example,  a circle is just a 
> segment with one end hooked to the other. Vertices are not duplicated, or 
> mapped to each other. This makes topology simple and easy to implement. 
> However, then geometry is more complicated. What Plex does is allow 
> coordinates to be represented by a discontinuous field taking values on 
> cells, in addition to vertices. In our circle example, each cells near the 
> cut will have 2 coordinates, one for each vertex, but they will not agree 
> across the cut. If you define a periodic domain, then Plex can construct this 
> coordinate field automatically using DMPlexLocalize(). These DG coordinates 
> are then used by the integration routines.

Ok, I think I understand the concept. DMPlex reads information about both 
topology and coordinates from the .msh file. Creating a periodic mesh in gmsh 
then should allow DMPlex to identify the periodic boundaries as the “cut” and 
build the mesh topology accordingly. Coordinate information is handled 
separately.
That means, as Lawrence suggested, building periodic meshes in gmsh and reading 
them in to petsc’s DMPlex should indeed “just work”.  (From the user 
perspective). The only extra step is to call DMLocalizeCoordinates() after 
DMPlexReadFromFile(). Sorry to reiterate, I am just trying to make sense of 
this. 
>  
> From my side, the issue is to do with mapping that coordinate field into one 
> that I understand (within Firedrake). You may not have this problem.
> 
> Firedrake uses its own coordinate mapping and integration routines, so they 
> must manage the second part independently. I hope to get change this slightly 
> soon by making the Firedrake representation a DMField, so that it looks the 
> same to Plex.
> 
>   Thanks,
> 
> Matt
>  
> Thanks,
> 
> Lawrence
> 
> 
> -- 
> 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
> 
> https://www.cse.buffalo.edu/~knepley/ 
> To read periodic meshes from GMSH, you need to use the option 
> -dm_plex_gmsh_periodic and DMPlexCreateFromFile

Ahh, thanks. I was missing the option " -dm_plex_gmsh_periodic “. But using 
this option I now generate a segmentation fault error when calling VecView() on 
the solution vector with vtk and hdf5 viewers. Any suggestions?

 
> 
> See  src/dm/impls/plex/examples/tests/ex1.c. An example runs
> 
> $ ./ex1 -filename 
> ${PETSC_DIR}/share/petsc/datafiles/meshes/cube_periodic_bin.msh 
> -dm_plex_gmsh_periodic -dm_view ::ascii_info_detail -interpolate -test_shape
> 
> generating periodic meshes in gmsh may be tricky, Lisandro for sure may 
> advice.

Thanks,
Max

Re: [petsc-users] Problems about Compiling Multifile Program

[Balay, Satish]

Please send complete logs. [i.e the complete compile output from your makefile. 
Also send your makefile]

If your c++ code uses .cxx - it should get compiled with a c++ compiler.

You might benefit by change CLINKER to CXXLINKER

Satish

On Wed, 24 Oct 2018, Yingjie Wu wrote:

> Thank you very much for your reply.
> 
> Because the files I added are written in c++, including definitions of
> classes and functions. I encountered a mistake when I used your method.
> myprogram: class1.o myprogram.o
>${CLINKER} -o myprogram class1.o myprogram.o ${PETSC_LIB}
> Error information is : unknown type name 'class'
> It seems that I did not compile using the c++ compiler.
> 
> How should I set it up so that it can compile successfully?
> 
> Thanks,
> Yingjie
> 
> Matthew Knepley  于2018年10月24日周三 上午1:46写道：
> 
> > On Tue, Oct 23, 2018 at 11:37 AM Yingjie Wu  wrote:
> >
> >> Dear Petsc developer:
> >> Hi,
> >> Thank you very much for your continuous help, I recently encountered some
> >> difficulties in developing programs on Petsc.
> >>
> >> 1. I want to use my class definition (class1. h) and class functions
> >> (class1. cpp) files in my Petsc program (myprogram. c) and compile my
> >> program. I have written other program before:
> >> g++ -c class1.cpp myprogram.c
> >> g++ -o myprogram *.o
> >> I also have some knowledge of Makefile, but I don't know how to modify it
> >> to achieve compilation. I really want your help.
> >>
> >
> > With PETSc makefiles, we have automatic rules, so you can use
> >
> > myprogram: class1.o myprogram.o
> >${CLINKER} -o myprogram class1.o myprogram.o ${PETSC_LIB}
> >
> >
> >> 2.  In my class definition and function implementation program, I want to
> >> use Petsc data structures, such as PetscInt, PetscReal. In order to
> >> maintain consistency of the main program. What kind of header files should
> >> I add to my code, or do I need to add any code?
> >>
> >
> > The easiest thing to do is just
> >
> > #include 
> >
> >   Thanks,
> >
> > Matt
> >
> >
> >> Thanks,
> >> Yingjie
> >>
> >
> >
> > --
> > 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
> >
> > https://www.cse.buffalo.edu/~knepley/
> > 
> >
> 

Re: [petsc-users] Problems about Compiling Multifile Program

[Yingjie Wu]

Thank you very much for your reply.

Because the files I added are written in c++, including definitions of
classes and functions. I encountered a mistake when I used your method.
myprogram: class1.o myprogram.o
   ${CLINKER} -o myprogram class1.o myprogram.o ${PETSC_LIB}
Error information is : unknown type name 'class'
It seems that I did not compile using the c++ compiler.

How should I set it up so that it can compile successfully?

Thanks,
Yingjie

Matthew Knepley  于2018年10月24日周三 上午1:46写道：

> On Tue, Oct 23, 2018 at 11:37 AM Yingjie Wu  wrote:
>
>> Dear Petsc developer:
>> Hi,
>> Thank you very much for your continuous help, I recently encountered some
>> difficulties in developing programs on Petsc.
>>
>> 1. I want to use my class definition (class1. h) and class functions
>> (class1. cpp) files in my Petsc program (myprogram. c) and compile my
>> program. I have written other program before:
>> g++ -c class1.cpp myprogram.c
>> g++ -o myprogram *.o
>> I also have some knowledge of Makefile, but I don't know how to modify it
>> to achieve compilation. I really want your help.
>>
>
> With PETSc makefiles, we have automatic rules, so you can use
>
> myprogram: class1.o myprogram.o
>${CLINKER} -o myprogram class1.o myprogram.o ${PETSC_LIB}
>
>
>> 2.  In my class definition and function implementation program, I want to
>> use Petsc data structures, such as PetscInt, PetscReal. In order to
>> maintain consistency of the main program. What kind of header files should
>> I add to my code, or do I need to add any code?
>>
>
> The easiest thing to do is just
>
> #include 
>
>   Thanks,
>
> Matt
>
>
>> Thanks,
>> Yingjie
>>
>
>
> --
> 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
>
> https://www.cse.buffalo.edu/~knepley/
> 
>

[petsc-users] Problems with a Laplacian problem with Neumann boundary condition

[SONG Pengyang]

Hi all,
I am trying to change PETSC KSP tutorial ex34 to fortran code ex34f and I 
referred to ex22f, but it somehow seems wrong in my code.
Important error information points to MatSetValuesStencil() [argument out of 
range].
Can anyone look over my code and give any advice?
Many thanks.

SONG Pengyang （宋朋洋）
MSc candidate in Physical Oceanography
College of Oceanic and Atmospheric Sciences
Ocean University of China
No. 238 Songling Road, Qingdao 266100, Shandong, China
TEL: (+86)15762264643



ex34f.F90
Description: ex34f.F90

Re: [petsc-users] PetscInt overflow

[Matthew Knepley]

On Wed, Oct 24, 2018 at 10:03 AM Jan Grießer 
wrote:

> This is the error message i get from my program:
> Shape of the dynamical matrix:  (150, 150)
>

Petsc installs a signal handler, so there should be a PETSc-specific
message before this one. Is something eating
your output?

  Thanks,

 Matt


>
> ===
> =   BAD TERMINATION OF ONE OF YOUR APPLICATION PROCESSES
> =   PID 122676 RUNNING AT n3512.nemo.privat
> =   EXIT CODE: 9
> =   CLEANING UP REMAINING PROCESSES
> =   YOU CAN IGNORE THE BELOW CLEANUP MESSAGES
>
> ===
>Intel(R) MPI Library troubleshooting guide:
>   https://software.intel.com/node/561764
>
> ===
>
>
> Am Mi., 24. Okt. 2018 um 16:01 Uhr schrieb Matthew Knepley <
> knep...@gmail.com>:
>
>> On Wed, Oct 24, 2018 at 9:38 AM Jan Grießer 
>> wrote:
>>
>>> Hey,
>>> i tried to run my program as you said with -bv_type vecs and/or -bv_type
>>> mat, but instead of the PETScInt overflow i now get an MPI Error 9
>>>
>>
>> Send the actual error.
>>
>>   Thanks,
>>
>> Matt
>>
>>
>>> message, which i assume (after googling a little bit around) should be a
>>> memory problem. I tried to run it also on slightly bigger compute nodes on
>>> our cluster with 20 cores with each 12 GB and 24 GB but the problem still
>>> remains. The actual limit appears to be a 1.5 Million x 1.5 Million where
>>> i  searched for NEV = 1500 and MPD = 500/ 200 eigenvalues.
>>> Do you have maybe an idea what the error could be? I appended also the
>>> python method i used to solve the problem. I also tried to solve the
>>> problem with spectrum solving but the error message remains the same.
>>>
>>> def solve_eigensystem(DynMatrix_nn, NEV, MPD, Dimension):
>>> # Create the solver
>>> # E is used as an acronym for the EPS solver (EPS = Eigenvalue problem
>>> solver)
>>> E = SLEPc.EPS().create()
>>>
>>> # Set the preconditioner
>>> pc = PETSc.PC().create()
>>> pc.setType(pc.Type.BJACOBI)
>>>
>>> # Set the linear solver
>>> # Create the KSP object
>>> ksp = PETSc.KSP().create()
>>> # Create the solver, in this case GMRES
>>> ksp.setType(ksp.Type.GMRES)
>>> # Set the tolerances of the GMRES solver
>>> # Link it to the PC
>>> ksp.setPC(pc)
>>>
>>> # Set up the spectral transformations
>>> st = SLEPc.ST().create()
>>> st.setType("shift")
>>> st.setKSP(ksp)
>>> # MPD stands for "maximum projected dimension". It has to due with
>>> computational cost, please read Chap. 2.6.5 of SLEPc docu for
>>> # an explanation. At the moment mpd is only a guess
>>> E.setDimensions(nev=NEV, mpd = MPD)
>>> # Eigenvalues should be real, therefore we start to order them from the
>>> smallest real value |l.real|
>>> E.setWhichEigenpairs(E.Which.SMALLEST_REAL)
>>> # Since the dynamical matrix is symmetric and real it is hermitian
>>> E.setProblemType(SLEPc.EPS.ProblemType.HEP)
>>> # Use the Krylov Schur method to solve the eigenvalue problem
>>> E.setType(E.Type.KRYLOVSCHUR)
>>> # Set the convergence criterion to relative to the eigenvalue and the
>>> maximal number of iterations
>>> E.setConvergenceTest(E.Conv.REL)
>>> E.setTolerances(tol = 1e-8, max_it = 5000)
>>> # Set the matrix in order to solve
>>> E.setOperators(DynMatrix_nn, None)
>>> # Sets EPS options from the options database. This routine must be
>>> called before `setUp()` if the user is to be allowed to set dthe solver
>>> type.
>>> E.setFromOptions()
>>> # Sets up all the internal data structures necessary for the execution
>>> of the eigensolver.
>>> E.setUp()
>>>
>>> # Solve eigenvalue problem
>>> E.solve()
>>>
>>> Print = PETSc.Sys.Print
>>>
>>> Print()
>>> Print("")
>>> Print("***SLEPc Solution Results***")
>>> Print("")
>>>
>>> its = E.getIterationNumber()
>>> Print("Number of iterations of the method: ", its)
>>> eps_type = E.getType()
>>> Print("Solution method: ", eps_type)
>>> nev, ncv, mpd = E.getDimensions()
>>> Print("Number of requested eigenvalues: ", nev)
>>> Print("Number of computeded eigenvectors: ", ncv)
>>> tol, maxit = E.getTolerances()
>>> Print("Stopping condition: (tol, maxit)", (tol, maxit))
>>> # Get the type of convergence
>>> conv_test = E.getConvergenceTest()
>>> Print("Selected convergence test: ", conv_test)
>>> # Get the used spectral transformation
>>> get_st = E.getST()
>>> Print("Selected spectral transformation: ", get_st)
>>> # Get the applied direct solver
>>> get_ksp = E.getDS()
>>> Print("Selected direct solver: ", get_ksp)
>>> nconv = E.getConverged()
>>> Print("Number of converged eigenpairs: ", nconv)
>>> .
>>>
>>>
>>>
>>> Am Fr., 19. Okt. 2018 um 21:00 Uhr schrieb Smith, Barry F. <
>>> bsm...@mcs.anl.gov>:
>>>


 > On Oct 19, 2018, at 7:56 AM, Zhang, Junchao 
 wrote:
 >
 >
 > On Fri, Oct 19, 2018 

Re: [petsc-users] PetscInt overflow

[Jose E. Roman]

The program seems correct, although in your case you don't need any KSP or 
preconditioner, so you could remove all code related to PC, KSP and ST - 
anyway, this should not affect memory consumption.

Add the option -eps_view_pre and send the output.
Jose


> El 24 oct 2018, a las 16:03, Jan Grießer  
> escribió:
> 
> This is the error message i get from my program: 
> Shape of the dynamical matrix:  (150, 150)
> 
> ===
> =   BAD TERMINATION OF ONE OF YOUR APPLICATION PROCESSES
> =   PID 122676 RUNNING AT n3512.nemo.privat
> =   EXIT CODE: 9
> =   CLEANING UP REMAINING PROCESSES
> =   YOU CAN IGNORE THE BELOW CLEANUP MESSAGES
> ===
>Intel(R) MPI Library troubleshooting guide:
>   https://software.intel.com/node/561764
> ===
> 
> 
> Am Mi., 24. Okt. 2018 um 16:01 Uhr schrieb Matthew Knepley 
> :
> On Wed, Oct 24, 2018 at 9:38 AM Jan Grießer  
> wrote:
> Hey,
> i tried to run my program as you said with -bv_type vecs and/or -bv_type mat, 
> but instead of the PETScInt overflow i now get an MPI Error 9
> 
> Send the actual error.
> 
>   Thanks,
> 
> Matt
>  
> message, which i assume (after googling a little bit around) should be a 
> memory problem. I tried to run it also on slightly bigger compute nodes on 
> our cluster with 20 cores with each 12 GB and 24 GB but the problem still 
> remains. The actual limit appears to be a 1.5 Million x 1.5 Million where i  
> searched for NEV = 1500 and MPD = 500/ 200 eigenvalues. 
> Do you have maybe an idea what the error could be? I appended also the python 
> method i used to solve the problem. I also tried to solve the problem with 
> spectrum solving but the error message remains the same. 
> 
> def solve_eigensystem(DynMatrix_nn, NEV, MPD, Dimension):
>   # Create the solver 
>   # E is used as an acronym for the EPS solver (EPS = Eigenvalue problem 
> solver)
>   E = SLEPc.EPS().create()
> 
>   # Set the preconditioner 
>   pc = PETSc.PC().create()
>   pc.setType(pc.Type.BJACOBI)
> 
>   # Set the linear solver 
>   # Create the KSP object
>   ksp = PETSc.KSP().create()
>   # Create the solver, in this case GMRES
>   ksp.setType(ksp.Type.GMRES)
>   # Set the tolerances of the GMRES solver
> # Link it to the PC
>   ksp.setPC(pc)
> 
>   # Set up the spectral transformations 
>   st = SLEPc.ST().create()
>   st.setType("shift")
>   st.setKSP(ksp)
>   
>   # MPD stands for "maximum projected dimension". It has to due with 
> computational cost, please read Chap. 2.6.5 of SLEPc docu for 
>   # an explanation. At the moment mpd is only a guess
>   E.setDimensions(nev=NEV, mpd = MPD)
>   # Eigenvalues should be real, therefore we start to order them from the 
> smallest real value |l.real|
>   E.setWhichEigenpairs(E.Which.SMALLEST_REAL)
>   # Since the dynamical matrix is symmetric and real it is hermitian 
>   E.setProblemType(SLEPc.EPS.ProblemType.HEP)
>   # Use the Krylov Schur method to solve the eigenvalue problem
>   E.setType(E.Type.KRYLOVSCHUR)
>   # Set the convergence criterion to relative to the eigenvalue and the 
> maximal number of iterations
>   E.setConvergenceTest(E.Conv.REL)
>   E.setTolerances(tol = 1e-8, max_it = 5000)
>   # Set the matrix in order to solve 
>   E.setOperators(DynMatrix_nn, None)
>   # Sets EPS options from the options database. This routine must be 
> called before `setUp()` if the user is to be allowed to set dthe solver type.
>   E.setFromOptions()
>   # Sets up all the internal data structures necessary for the execution 
> of the eigensolver.
>   E.setUp()
> 
>   # Solve eigenvalue problem 
>   E.solve()
> 
>   Print = PETSc.Sys.Print
> 
>   Print()
>   Print("")
>   Print("***SLEPc Solution Results***")
>   Print("")
> 
>   its = E.getIterationNumber()
>   Print("Number of iterations of the method: ", its)
>   eps_type = E.getType()
>   Print("Solution method: ", eps_type)
>   nev, ncv, mpd = E.getDimensions()
>   Print("Number of requested eigenvalues: ", nev)
>   Print("Number of computeded eigenvectors: ", ncv)
>   tol, maxit = E.getTolerances()
>   Print("Stopping condition: (tol, maxit)", (tol, maxit))
>   # Get the type of convergence
>   conv_test = E.getConvergenceTest()
>   Print("Selected convergence test: ", conv_test)
>   # Get the used spectral transformation
>   get_st = E.getST()
>   Print("Selected spectral transformation: ", get_st)
>   # Get the applied direct solver
>   get_ksp = E.getDS()
>   Print("Selected direct solver: ", get_ksp)
>   

Re: [petsc-users] PetscInt overflow

[Jan Grießer]

This is the error message i get from my program:
Shape of the dynamical matrix:  (150, 150)

===
=   BAD TERMINATION OF ONE OF YOUR APPLICATION PROCESSES
=   PID 122676 RUNNING AT n3512.nemo.privat
=   EXIT CODE: 9
=   CLEANING UP REMAINING PROCESSES
=   YOU CAN IGNORE THE BELOW CLEANUP MESSAGES
===
   Intel(R) MPI Library troubleshooting guide:
  https://software.intel.com/node/561764
===


Am Mi., 24. Okt. 2018 um 16:01 Uhr schrieb Matthew Knepley <
knep...@gmail.com>:

> On Wed, Oct 24, 2018 at 9:38 AM Jan Grießer 
> wrote:
>
>> Hey,
>> i tried to run my program as you said with -bv_type vecs and/or -bv_type
>> mat, but instead of the PETScInt overflow i now get an MPI Error 9
>>
>
> Send the actual error.
>
>   Thanks,
>
> Matt
>
>
>> message, which i assume (after googling a little bit around) should be a
>> memory problem. I tried to run it also on slightly bigger compute nodes on
>> our cluster with 20 cores with each 12 GB and 24 GB but the problem still
>> remains. The actual limit appears to be a 1.5 Million x 1.5 Million where
>> i  searched for NEV = 1500 and MPD = 500/ 200 eigenvalues.
>> Do you have maybe an idea what the error could be? I appended also the
>> python method i used to solve the problem. I also tried to solve the
>> problem with spectrum solving but the error message remains the same.
>>
>> def solve_eigensystem(DynMatrix_nn, NEV, MPD, Dimension):
>> # Create the solver
>> # E is used as an acronym for the EPS solver (EPS = Eigenvalue problem
>> solver)
>> E = SLEPc.EPS().create()
>>
>> # Set the preconditioner
>> pc = PETSc.PC().create()
>> pc.setType(pc.Type.BJACOBI)
>>
>> # Set the linear solver
>> # Create the KSP object
>> ksp = PETSc.KSP().create()
>> # Create the solver, in this case GMRES
>> ksp.setType(ksp.Type.GMRES)
>> # Set the tolerances of the GMRES solver
>> # Link it to the PC
>> ksp.setPC(pc)
>>
>> # Set up the spectral transformations
>> st = SLEPc.ST().create()
>> st.setType("shift")
>> st.setKSP(ksp)
>> # MPD stands for "maximum projected dimension". It has to due with
>> computational cost, please read Chap. 2.6.5 of SLEPc docu for
>> # an explanation. At the moment mpd is only a guess
>> E.setDimensions(nev=NEV, mpd = MPD)
>> # Eigenvalues should be real, therefore we start to order them from the
>> smallest real value |l.real|
>> E.setWhichEigenpairs(E.Which.SMALLEST_REAL)
>> # Since the dynamical matrix is symmetric and real it is hermitian
>> E.setProblemType(SLEPc.EPS.ProblemType.HEP)
>> # Use the Krylov Schur method to solve the eigenvalue problem
>> E.setType(E.Type.KRYLOVSCHUR)
>> # Set the convergence criterion to relative to the eigenvalue and the
>> maximal number of iterations
>> E.setConvergenceTest(E.Conv.REL)
>> E.setTolerances(tol = 1e-8, max_it = 5000)
>> # Set the matrix in order to solve
>> E.setOperators(DynMatrix_nn, None)
>> # Sets EPS options from the options database. This routine must be called
>> before `setUp()` if the user is to be allowed to set dthe solver type.
>> E.setFromOptions()
>> # Sets up all the internal data structures necessary for the execution of
>> the eigensolver.
>> E.setUp()
>>
>> # Solve eigenvalue problem
>> E.solve()
>>
>> Print = PETSc.Sys.Print
>>
>> Print()
>> Print("")
>> Print("***SLEPc Solution Results***")
>> Print("")
>>
>> its = E.getIterationNumber()
>> Print("Number of iterations of the method: ", its)
>> eps_type = E.getType()
>> Print("Solution method: ", eps_type)
>> nev, ncv, mpd = E.getDimensions()
>> Print("Number of requested eigenvalues: ", nev)
>> Print("Number of computeded eigenvectors: ", ncv)
>> tol, maxit = E.getTolerances()
>> Print("Stopping condition: (tol, maxit)", (tol, maxit))
>> # Get the type of convergence
>> conv_test = E.getConvergenceTest()
>> Print("Selected convergence test: ", conv_test)
>> # Get the used spectral transformation
>> get_st = E.getST()
>> Print("Selected spectral transformation: ", get_st)
>> # Get the applied direct solver
>> get_ksp = E.getDS()
>> Print("Selected direct solver: ", get_ksp)
>> nconv = E.getConverged()
>> Print("Number of converged eigenpairs: ", nconv)
>> .
>>
>>
>>
>> Am Fr., 19. Okt. 2018 um 21:00 Uhr schrieb Smith, Barry F. <
>> bsm...@mcs.anl.gov>:
>>
>>>
>>>
>>> > On Oct 19, 2018, at 7:56 AM, Zhang, Junchao 
>>> wrote:
>>> >
>>> >
>>> > On Fri, Oct 19, 2018 at 4:02 AM Jan Grießer <
>>> griesser@googlemail.com> wrote:
>>> > With more than 1 MPI process you mean i should use spectrum slicing in
>>> divide the full problem in smaller subproblems?
>>> > The --with-64-bit-indices is not a possibility for me since i
>>> configured petsc with mumps, which does not allow to use the 64-bit version
>>> (At least this 

Re: [petsc-users] PetscInt overflow

[Matthew Knepley]

On Wed, Oct 24, 2018 at 9:38 AM Jan Grießer 
wrote:

> Hey,
> i tried to run my program as you said with -bv_type vecs and/or -bv_type
> mat, but instead of the PETScInt overflow i now get an MPI Error 9
>

Send the actual error.

  Thanks,

Matt


> message, which i assume (after googling a little bit around) should be a
> memory problem. I tried to run it also on slightly bigger compute nodes on
> our cluster with 20 cores with each 12 GB and 24 GB but the problem still
> remains. The actual limit appears to be a 1.5 Million x 1.5 Million where
> i  searched for NEV = 1500 and MPD = 500/ 200 eigenvalues.
> Do you have maybe an idea what the error could be? I appended also the
> python method i used to solve the problem. I also tried to solve the
> problem with spectrum solving but the error message remains the same.
>
> def solve_eigensystem(DynMatrix_nn, NEV, MPD, Dimension):
> # Create the solver
> # E is used as an acronym for the EPS solver (EPS = Eigenvalue problem
> solver)
> E = SLEPc.EPS().create()
>
> # Set the preconditioner
> pc = PETSc.PC().create()
> pc.setType(pc.Type.BJACOBI)
>
> # Set the linear solver
> # Create the KSP object
> ksp = PETSc.KSP().create()
> # Create the solver, in this case GMRES
> ksp.setType(ksp.Type.GMRES)
> # Set the tolerances of the GMRES solver
> # Link it to the PC
> ksp.setPC(pc)
>
> # Set up the spectral transformations
> st = SLEPc.ST().create()
> st.setType("shift")
> st.setKSP(ksp)
> # MPD stands for "maximum projected dimension". It has to due with
> computational cost, please read Chap. 2.6.5 of SLEPc docu for
> # an explanation. At the moment mpd is only a guess
> E.setDimensions(nev=NEV, mpd = MPD)
> # Eigenvalues should be real, therefore we start to order them from the
> smallest real value |l.real|
> E.setWhichEigenpairs(E.Which.SMALLEST_REAL)
> # Since the dynamical matrix is symmetric and real it is hermitian
> E.setProblemType(SLEPc.EPS.ProblemType.HEP)
> # Use the Krylov Schur method to solve the eigenvalue problem
> E.setType(E.Type.KRYLOVSCHUR)
> # Set the convergence criterion to relative to the eigenvalue and the
> maximal number of iterations
> E.setConvergenceTest(E.Conv.REL)
> E.setTolerances(tol = 1e-8, max_it = 5000)
> # Set the matrix in order to solve
> E.setOperators(DynMatrix_nn, None)
> # Sets EPS options from the options database. This routine must be called
> before `setUp()` if the user is to be allowed to set dthe solver type.
> E.setFromOptions()
> # Sets up all the internal data structures necessary for the execution of
> the eigensolver.
> E.setUp()
>
> # Solve eigenvalue problem
> E.solve()
>
> Print = PETSc.Sys.Print
>
> Print()
> Print("")
> Print("***SLEPc Solution Results***")
> Print("")
>
> its = E.getIterationNumber()
> Print("Number of iterations of the method: ", its)
> eps_type = E.getType()
> Print("Solution method: ", eps_type)
> nev, ncv, mpd = E.getDimensions()
> Print("Number of requested eigenvalues: ", nev)
> Print("Number of computeded eigenvectors: ", ncv)
> tol, maxit = E.getTolerances()
> Print("Stopping condition: (tol, maxit)", (tol, maxit))
> # Get the type of convergence
> conv_test = E.getConvergenceTest()
> Print("Selected convergence test: ", conv_test)
> # Get the used spectral transformation
> get_st = E.getST()
> Print("Selected spectral transformation: ", get_st)
> # Get the applied direct solver
> get_ksp = E.getDS()
> Print("Selected direct solver: ", get_ksp)
> nconv = E.getConverged()
> Print("Number of converged eigenpairs: ", nconv)
> .
>
>
>
> Am Fr., 19. Okt. 2018 um 21:00 Uhr schrieb Smith, Barry F. <
> bsm...@mcs.anl.gov>:
>
>>
>>
>> > On Oct 19, 2018, at 7:56 AM, Zhang, Junchao 
>> wrote:
>> >
>> >
>> > On Fri, Oct 19, 2018 at 4:02 AM Jan Grießer <
>> griesser@googlemail.com> wrote:
>> > With more than 1 MPI process you mean i should use spectrum slicing in
>> divide the full problem in smaller subproblems?
>> > The --with-64-bit-indices is not a possibility for me since i
>> configured petsc with mumps, which does not allow to use the 64-bit version
>> (At least this was the error message when i tried to configure PETSc )
>> >
>> > MUMPS 5.1.2 manual chapter 2.4.2 says it supports "Selective 64-bit
>> integer feature" and "full 64-bit integer version" as well.
>>
>> They use to achieve this by compiling with special Fortran flags to
>> promote integers to 64 bit; this is too fragile for our taste so we never
>> hooked PETSc up wit it. If they have a dependable way of using 64 bit
>> integers we should add that to our mumps.py and test it.
>>
>>Barry
>>
>> >
>> > Am Mi., 17. Okt. 2018 um 18:24 Uhr schrieb Jose E. Roman <
>> jro...@dsic.upv.es>:
>> > To use BVVECS just add the command-line option -bv_type vecs
>> > This causes to use a separate Vec for each column, instead of a single
>> long Vec of size n*m. But it is considerably slower than the default.
>> >
>> > Anyway, for such large problems you 

Re: [petsc-users] [SLEPc] Number of iterations changes with MPI processes in Lanczos

[Jose E. Roman]

Everything seems correct. I don't know, maybe your problem is very sensitive? 
Is the eigenvalue tiny?
I would still try with Krylov-Schur.
Jose


> El 24 oct 2018, a las 14:59, Ale Foggia  escribió:
> 
> The functions called to set the solver are (in this order): EPSCreate(); 
> EPSSetOperators(); EPSSetProblemType(EPS_HEP); EPSSetType(EPSLANCZOS); 
> EPSSetWhichEigenpairs(EPS_SMALLEST_REAL); EPSSetFromOptions();
> 
> The output of -eps_view for each run is:
> =
> EPS Object: 960 MPI processes
>   type: lanczos
> LOCAL reorthogonalization
>   problem type: symmetric eigenvalue problem
>   selected portion of the spectrum: smallest real parts
>   number of eigenvalues (nev): 1
>   number of column vectors (ncv): 16
>   maximum dimension of projected problem (mpd): 16
>   maximum number of iterations: 291700777
>   tolerance: 1e-08
>   convergence test: relative to the eigenvalue
> BV Object: 960 MPI processes
>   type: svec
>   17 columns of global length 2333606220
>   vector orthogonalization method: modified Gram-Schmidt
>   orthogonalization refinement: if needed (eta: 0.7071)
>   block orthogonalization method: GS
>   doing matmult as a single matrix-matrix product
>   generating random vectors independent of the number of processes
> DS Object: 960 MPI processes
>   type: hep
>   parallel operation mode: REDUNDANT
>   solving the problem with: Implicit QR method (_steqr)
> ST Object: 960 MPI processes
>   type: shift
>   shift: 0.
>   number of matrices: 1
> =
> EPS Object: 1024 MPI processes
>   type: lanczos
> LOCAL reorthogonalization
>   problem type: symmetric eigenvalue problem
>   selected portion of the spectrum: smallest real parts
>   number of eigenvalues (nev): 1
>   number of column vectors (ncv): 16
>   maximum dimension of projected problem (mpd): 16
>   maximum number of iterations: 291700777
>   tolerance: 1e-08
>   convergence test: relative to the eigenvalue
> BV Object: 1024 MPI processes
>   type: svec
>   17 columns of global length 2333606220
>   vector orthogonalization method: modified Gram-Schmidt
>   orthogonalization refinement: if needed (eta: 0.7071)
>   block orthogonalization method: GS
>   doing matmult as a single matrix-matrix product
>   generating random vectors independent of the number of processes
> DS Object: 1024 MPI processes
>   type: hep
>   parallel operation mode: REDUNDANT
>   solving the problem with: Implicit QR method (_steqr)
> ST Object: 1024 MPI processes
>   type: shift
>   shift: 0.
>   number of matrices: 1
> =
> 
> I run again the same configurations and I got the same result in term of the 
> number of iterations.
> 
> I also tried the full reorthogonalization (always with the 
> -bv_reproducible_random option) but I still get different number of 
> iterations: for 960 procs I get 172 iters, and for 1024 I get 362 iters. The 
> -esp_view output for this case (only for 960 procs, the other one has the 
> same information -except the number of processes-) is:
> =
> EPS Object: 960 MPI processes
>   type: lanczos
> FULL reorthogonalization
>   problem type: symmetric eigenvalue problem
>   selected portion of the spectrum: smallest real parts
>   number of eigenvalues (nev): 1
>   number of column vectors (ncv): 16
>   maximum dimension of projected problem (mpd): 16
>   maximum number of iterations: 291700777
>   tolerance: 1e-08
>   convergence test: relative to the eigenvalue
> BV Object: 960 MPI processes
>   type: svec
>   17 columns of global length 2333606220
>   vector orthogonalization method: classical Gram-Schmidt
>   orthogonalization refinement: if needed (eta: 0.7071)
>   block orthogonalization method: GS
>   doing matmult as a single matrix-matrix product
>   generating random vectors independent of the number of processes
> DS Object: 960 MPI processes
>   type: hep
>   parallel operation mode: REDUNDANT
>   solving the problem with: Implicit QR method (_steqr)
> ST Object: 960 MPI processes
>   type: shift
>   shift: 0.
>   number of matrices: 1
> =
> 
> El mié., 24 oct. 2018 a las 10:52, Jose E. Roman () 
> escribió:
> This is very strange. Make sure you call EPSSetFromOptions in the code. Do 
> iteration counts change also for two different runs with the same number of 
> processes?
> Maybe Lanczos with default options is too sensitive (by default it does not 
> reorthogonalize). Suggest using Krylov-Schur or Lanczos with full 
> reorthogonalization (EPSLanczosSetReorthog).
> Also, send the output of -eps_view to see if there is anything abnormal.
> 
> Jose
> 
> 
> > El 24 oct 2018, a las 9:09, Ale Foggia  escribió:
> > 
> > I've tried the option that you gave me but I still get different number of 

Re: [petsc-users] PetscInt overflow

[Jan Grießer]

Hey,
i tried to run my program as you said with -bv_type vecs and/or -bv_type
mat, but instead of the PETScInt overflow i now get an MPI Error 9 message,
which i assume (after googling a little bit around) should be a memory
problem. I tried to run it also on slightly bigger compute nodes on our
cluster with 20 cores with each 12 GB and 24 GB but the problem still
remains. The actual limit appears to be a 1.5 Million x 1.5 Million where
i  searched for NEV = 1500 and MPD = 500/ 200 eigenvalues.
Do you have maybe an idea what the error could be? I appended also the
python method i used to solve the problem. I also tried to solve the
problem with spectrum solving but the error message remains the same.

def solve_eigensystem(DynMatrix_nn, NEV, MPD, Dimension):
# Create the solver
# E is used as an acronym for the EPS solver (EPS = Eigenvalue problem
solver)
E = SLEPc.EPS().create()

# Set the preconditioner
pc = PETSc.PC().create()
pc.setType(pc.Type.BJACOBI)

# Set the linear solver
# Create the KSP object
ksp = PETSc.KSP().create()
# Create the solver, in this case GMRES
ksp.setType(ksp.Type.GMRES)
# Set the tolerances of the GMRES solver
# Link it to the PC
ksp.setPC(pc)

# Set up the spectral transformations
st = SLEPc.ST().create()
st.setType("shift")
st.setKSP(ksp)
# MPD stands for "maximum projected dimension". It has to due with
computational cost, please read Chap. 2.6.5 of SLEPc docu for
# an explanation. At the moment mpd is only a guess
E.setDimensions(nev=NEV, mpd = MPD)
# Eigenvalues should be real, therefore we start to order them from the
smallest real value |l.real|
E.setWhichEigenpairs(E.Which.SMALLEST_REAL)
# Since the dynamical matrix is symmetric and real it is hermitian
E.setProblemType(SLEPc.EPS.ProblemType.HEP)
# Use the Krylov Schur method to solve the eigenvalue problem
E.setType(E.Type.KRYLOVSCHUR)
# Set the convergence criterion to relative to the eigenvalue and the
maximal number of iterations
E.setConvergenceTest(E.Conv.REL)
E.setTolerances(tol = 1e-8, max_it = 5000)
# Set the matrix in order to solve
E.setOperators(DynMatrix_nn, None)
# Sets EPS options from the options database. This routine must be called
before `setUp()` if the user is to be allowed to set dthe solver type.
E.setFromOptions()
# Sets up all the internal data structures necessary for the execution of
the eigensolver.
E.setUp()

# Solve eigenvalue problem
E.solve()

Print = PETSc.Sys.Print

Print()
Print("")
Print("***SLEPc Solution Results***")
Print("")

its = E.getIterationNumber()
Print("Number of iterations of the method: ", its)
eps_type = E.getType()
Print("Solution method: ", eps_type)
nev, ncv, mpd = E.getDimensions()
Print("Number of requested eigenvalues: ", nev)
Print("Number of computeded eigenvectors: ", ncv)
tol, maxit = E.getTolerances()
Print("Stopping condition: (tol, maxit)", (tol, maxit))
# Get the type of convergence
conv_test = E.getConvergenceTest()
Print("Selected convergence test: ", conv_test)
# Get the used spectral transformation
get_st = E.getST()
Print("Selected spectral transformation: ", get_st)
# Get the applied direct solver
get_ksp = E.getDS()
Print("Selected direct solver: ", get_ksp)
nconv = E.getConverged()
Print("Number of converged eigenpairs: ", nconv)
.



Am Fr., 19. Okt. 2018 um 21:00 Uhr schrieb Smith, Barry F. <
bsm...@mcs.anl.gov>:

>
>
> > On Oct 19, 2018, at 7:56 AM, Zhang, Junchao  wrote:
> >
> >
> > On Fri, Oct 19, 2018 at 4:02 AM Jan Grießer 
> wrote:
> > With more than 1 MPI process you mean i should use spectrum slicing in
> divide the full problem in smaller subproblems?
> > The --with-64-bit-indices is not a possibility for me since i configured
> petsc with mumps, which does not allow to use the 64-bit version (At least
> this was the error message when i tried to configure PETSc )
> >
> > MUMPS 5.1.2 manual chapter 2.4.2 says it supports "Selective 64-bit
> integer feature" and "full 64-bit integer version" as well.
>
> They use to achieve this by compiling with special Fortran flags to
> promote integers to 64 bit; this is too fragile for our taste so we never
> hooked PETSc up wit it. If they have a dependable way of using 64 bit
> integers we should add that to our mumps.py and test it.
>
>Barry
>
> >
> > Am Mi., 17. Okt. 2018 um 18:24 Uhr schrieb Jose E. Roman <
> jro...@dsic.upv.es>:
> > To use BVVECS just add the command-line option -bv_type vecs
> > This causes to use a separate Vec for each column, instead of a single
> long Vec of size n*m. But it is considerably slower than the default.
> >
> > Anyway, for such large problems you should consider using more than 1
> MPI process. In that case the error may disappear because the local size is
> smaller than 768000.
> >
> > Jose
> >
> >
> > > El 17 oct 2018, a las 17:58, Matthew Knepley 
> escribió:
> > >
> > > On Wed, Oct 17, 2018 at 11:54 AM Jan Grießer <
> griesser@googlemail.com> wrote:
> > > Hi all,
> > > i am 

Re: [petsc-users] Periodic BC in petsc DMPlex / firedrake

[Stefano Zampini]

To read periodic meshes from GMSH, you need to use the option
-dm_plex_gmsh_periodic and DMPlexCreateFromFile
See  src/dm/impls/plex/examples/tests/ex1.c. An example runs

$ ./ex1 -filename
${PETSC_DIR}/share/petsc/datafiles/meshes/cube_periodic_bin.msh
-dm_plex_gmsh_periodic -dm_view ::ascii_info_detail -interpolate -test_shape

generating periodic meshes in gmsh may be tricky, Lisandro for sure may
advice.

Il giorno mer 24 ott 2018 alle ore 13:51 Matthew Knepley 
ha scritto:

> On Wed, Oct 24, 2018 at 6:29 AM Lawrence Mitchell  wrote:
>
>> Hi Max,
>>
>> (I'm cc'ing in the petsc-users mailing list which may have more advice,
>> if you are using PETSc you should definitely subscribe!
>>
>> > On 24 Oct 2018, at 09:27, Maximilian Hartig 
>> wrote:
>> >
>> > Hello Lawrence,
>> >
>> > sorry to message you out of the blue. My name is Max and I found your
>> post on GitHub (https://github.com/firedrakeproject/firedrake/issues/1246
>> ) on DMPlex being able to read periodic gmsh files. I am currently trying
>> to do just that (creating a periodic DMPlex mesh with gmsh) in the context
>> of my PhD work. So far I haven’ t found any documentation on the periodic
>> BC’s with DMPlex and gmsh in the official petsc documentation.
>> > I was wondering whether you’d be so kind as to point me in a general
>> direction concerning how to achieve this. You seem experienced in using
>> petsc and I would greatly appreciate your help.
>>
>>
>> I think the answer is "it depends". If you're just using DMPlex directly
>> and all the of the functionality with PetscDS, then I /think/ that reading
>> periodic meshes via gmsh (assuming you're using the appropriate gmsh mesh
>> format [v2]) "just works".
>>
>
> There are two phases here: topological and geometric. DMPlex represents
> the periodic topological entity directly. For example,  a circle is just a
> segment with one end hooked to the other. Vertices are not duplicated, or
> mapped to each other. This makes topology simple and easy to implement.
> However, then geometry is more complicated. What Plex does is allow
> coordinates to be represented by a discontinuous field taking values on
> cells, in addition to vertices. In our circle example, each cells near the
> cut will have 2 coordinates, one for each vertex, but they will not agree
> across the cut. If you define a periodic domain, then Plex can construct
> this coordinate field automatically using DMPlexLocalize(). These DG
> coordinates are then used by the integration routines.
>
>
>> From my side, the issue is to do with mapping that coordinate field into
>> one that I understand (within Firedrake). You may not have this problem.
>>
>
> Firedrake uses its own coordinate mapping and integration routines, so
> they must manage the second part independently. I hope to get change this
> slightly soon by making the Firedrake representation a DMField, so that it
> looks the same to Plex.
>
>   Thanks,
>
> Matt
>
>
>> Thanks,
>>
>> Lawrence
>
>
>
> --
> 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
>
> https://www.cse.buffalo.edu/~knepley/
> 
>


-- 
Stefano

Re: [petsc-users] [SLEPc] Number of iterations changes with MPI processes in Lanczos

[Ale Foggia]

The functions called to set the solver are (in this order): EPSCreate();
EPSSetOperators(); EPSSetProblemType(EPS_HEP); EPSSetType(EPSLANCZOS);
EPSSetWhichEigenpairs(EPS_SMALLEST_REAL); EPSSetFromOptions();

The output of -eps_view for each run is:
=
EPS Object: 960 MPI processes
  type: lanczos
LOCAL reorthogonalization
  problem type: symmetric eigenvalue problem
  selected portion of the spectrum: smallest real parts
  number of eigenvalues (nev): 1
  number of column vectors (ncv): 16
  maximum dimension of projected problem (mpd): 16
  maximum number of iterations: 291700777
  tolerance: 1e-08
  convergence test: relative to the eigenvalue
BV Object: 960 MPI processes
  type: svec
  17 columns of global length 2333606220
  vector orthogonalization method: modified Gram-Schmidt
  orthogonalization refinement: if needed (eta: 0.7071)
  block orthogonalization method: GS
  doing matmult as a single matrix-matrix product
  generating random vectors independent of the number of processes
DS Object: 960 MPI processes
  type: hep
  parallel operation mode: REDUNDANT
  solving the problem with: Implicit QR method (_steqr)
ST Object: 960 MPI processes
  type: shift
  shift: 0.
  number of matrices: 1
=
EPS Object: 1024 MPI processes
  type: lanczos
LOCAL reorthogonalization
  problem type: symmetric eigenvalue problem
  selected portion of the spectrum: smallest real parts
  number of eigenvalues (nev): 1
  number of column vectors (ncv): 16
  maximum dimension of projected problem (mpd): 16
  maximum number of iterations: 291700777
  tolerance: 1e-08
  convergence test: relative to the eigenvalue
BV Object: 1024 MPI processes
  type: svec
  17 columns of global length 2333606220
  vector orthogonalization method: modified Gram-Schmidt
  orthogonalization refinement: if needed (eta: 0.7071)
  block orthogonalization method: GS
  doing matmult as a single matrix-matrix product
  generating random vectors independent of the number of processes
DS Object: 1024 MPI processes
  type: hep
  parallel operation mode: REDUNDANT
  solving the problem with: Implicit QR method (_steqr)
ST Object: 1024 MPI processes
  type: shift
  shift: 0.
  number of matrices: 1
=

I run again the same configurations and I got the same result in term of
the number of iterations.

I also tried the full reorthogonalization (always with the
-bv_reproducible_random option) but I still get different number of
iterations: for 960 procs I get 172 iters, and for 1024 I get 362 iters.
The -esp_view output for this case (only for 960 procs, the other one has
the same information -except the number of processes-) is:
=
EPS Object: 960 MPI processes
  type: lanczos
FULL reorthogonalization
  problem type: symmetric eigenvalue problem
  selected portion of the spectrum: smallest real parts
  number of eigenvalues (nev): 1
  number of column vectors (ncv): 16
  maximum dimension of projected problem (mpd): 16
  maximum number of iterations: 291700777
  tolerance: 1e-08
  convergence test: relative to the eigenvalue
BV Object: 960 MPI processes
  type: svec
  17 columns of global length 2333606220
  vector orthogonalization method: classical Gram-Schmidt
  orthogonalization refinement: if needed (eta: 0.7071)
  block orthogonalization method: GS
  doing matmult as a single matrix-matrix product
  generating random vectors independent of the number of processes
DS Object: 960 MPI processes
  type: hep
  parallel operation mode: REDUNDANT
  solving the problem with: Implicit QR method (_steqr)
ST Object: 960 MPI processes
  type: shift
  shift: 0.
  number of matrices: 1
=

El mié., 24 oct. 2018 a las 10:52, Jose E. Roman ()
escribió:

> This is very strange. Make sure you call EPSSetFromOptions in the code. Do
> iteration counts change also for two different runs with the same number of
> processes?
> Maybe Lanczos with default options is too sensitive (by default it does
> not reorthogonalize). Suggest using Krylov-Schur or Lanczos with full
> reorthogonalization (EPSLanczosSetReorthog).
> Also, send the output of -eps_view to see if there is anything abnormal.
>
> Jose
>
>
> > El 24 oct 2018, a las 9:09, Ale Foggia  escribió:
> >
> > I've tried the option that you gave me but I still get different number
> of iterations when changing the number of MPI processes: I did 960 procs
> and 1024 procs and I got 435 and 176 iterations, respectively.
> >
> > El mar., 23 oct. 2018 a las 16:48, Jose E. Roman ()
> escribió:
> >
> >
> > > El 23 oct 2018, a las 15:46, Ale Foggia  escribió:
> > >
> > >
> > >
> > > El mar., 23 oct. 2018 a las 15:33, Jose E. Roman ()
> escribió:
> > >
> > >
> > > > El 23 oct 2018, a las 15:17, Ale Foggia 
> 

Re: [petsc-users] Periodic BC in petsc DMPlex / firedrake

[Matthew Knepley]

On Wed, Oct 24, 2018 at 6:29 AM Lawrence Mitchell  wrote:

> Hi Max,
>
> (I'm cc'ing in the petsc-users mailing list which may have more advice, if
> you are using PETSc you should definitely subscribe!
>
> > On 24 Oct 2018, at 09:27, Maximilian Hartig 
> wrote:
> >
> > Hello Lawrence,
> >
> > sorry to message you out of the blue. My name is Max and I found your
> post on GitHub (https://github.com/firedrakeproject/firedrake/issues/1246
> ) on DMPlex being able to read periodic gmsh files. I am currently trying
> to do just that (creating a periodic DMPlex mesh with gmsh) in the context
> of my PhD work. So far I haven’ t found any documentation on the periodic
> BC’s with DMPlex and gmsh in the official petsc documentation.
> > I was wondering whether you’d be so kind as to point me in a general
> direction concerning how to achieve this. You seem experienced in using
> petsc and I would greatly appreciate your help.
>
>
> I think the answer is "it depends". If you're just using DMPlex directly
> and all the of the functionality with PetscDS, then I /think/ that reading
> periodic meshes via gmsh (assuming you're using the appropriate gmsh mesh
> format [v2]) "just works".
>

There are two phases here: topological and geometric. DMPlex represents the
periodic topological entity directly. For example,  a circle is just a
segment with one end hooked to the other. Vertices are not duplicated, or
mapped to each other. This makes topology simple and easy to implement.
However, then geometry is more complicated. What Plex does is allow
coordinates to be represented by a discontinuous field taking values on
cells, in addition to vertices. In our circle example, each cells near the
cut will have 2 coordinates, one for each vertex, but they will not agree
across the cut. If you define a periodic domain, then Plex can construct
this coordinate field automatically using DMPlexLocalize(). These DG
coordinates are then used by the integration routines.


> From my side, the issue is to do with mapping that coordinate field into
> one that I understand (within Firedrake). You may not have this problem.
>

Firedrake uses its own coordinate mapping and integration routines, so they
must manage the second part independently. I hope to get change this
slightly soon by making the Firedrake representation a DMField, so that it
looks the same to Plex.

  Thanks,

Matt


> Thanks,
>
> Lawrence



-- 
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

https://www.cse.buffalo.edu/~knepley/ 

Re: [petsc-users] Periodic BC in petsc DMPlex / firedrake

[Lawrence Mitchell]

Hi Max,

(I'm cc'ing in the petsc-users mailing list which may have more advice, if you 
are using PETSc you should definitely subscribe!

> On 24 Oct 2018, at 09:27, Maximilian Hartig  wrote:
> 
> Hello Lawrence,
> 
> sorry to message you out of the blue. My name is Max and I found your post on 
> GitHub (https://github.com/firedrakeproject/firedrake/issues/1246 ) on DMPlex 
> being able to read periodic gmsh files. I am currently trying to do just that 
> (creating a periodic DMPlex mesh with gmsh) in the context of my PhD work. So 
> far I haven’ t found any documentation on the periodic BC’s with DMPlex and 
> gmsh in the official petsc documentation. 
> I was wondering whether you’d be so kind as to point me in a general 
> direction concerning how to achieve this. You seem experienced in using petsc 
> and I would greatly appreciate your help. 


I think the answer is "it depends". If you're just using DMPlex directly and 
all the of the functionality with PetscDS, then I /think/ that reading periodic 
meshes via gmsh (assuming you're using the appropriate gmsh mesh format [v2]) 
"just works".

From my side, the issue is to do with mapping that coordinate field into one 
that I understand (within Firedrake). You may not have this problem.

Thanks,

Lawrence

Re: [petsc-users] Using BDDC preconditioner for assembled matrices

[Stefano Zampini]

Abdullah,

The "Neumann" problems Jed is referring to result from assembling your
problem on each subdomain ( = MPI process) separately.
Assuming you are using FEM, these problems have been historically  named
"Neumann" as they correspond to a problem with natural boundary conditions
(Neumann bc for Poisson).
Note that in PETSc the subdomain decomposition is associated with the mesh
decomposition.

When converting from an assembled AIJ matrix to a MATIS format, such
"Neumann" information is lost.
You can disassemble an AIJ matrix, in the sense that you can find local
matrices A_j such that A = \sum_j R^T_j A_j R_j (as it is done in ex72.c),
but you cannot guarantee (unless if you solve an optimization problem) that
the disassembling will produce subdomain Neumann problems that are
consistent with your FEM problem.

I have added such disassembling code a few months ago, just to have another
alternative for preconditioning AIJ matrices in PETSc; there are few tweaks
one can do to improve the quality of the disassembling, but I discourage
its usage unless you don't have access to the FEM assembly code.

With that said, what problem are you trying to solve? Are you using DMDA or
DMPlex? What are the results you obtained with using the automatic
disassembling?

Il giorno mer 24 ott 2018 alle ore 08:14 Abdullah Ali Sivas <
abdullahasi...@gmail.com> ha scritto:

> Hi Jed,
>
> Thanks for your reply. The assembled matrix I have corresponds to the full
> problem on the full mesh. There are no "Neumann" problems (or any sort of
> domain decomposition) defined in the code generates the matrix. However, I
> think assembling the full problem is equivalent to implicitly assembling
> the "Neumann" problems, since the system can be partitioned as;
>
> [A_{LL} | A_{LI}]  [u_L] [F]
> ---|  = -
> [A_{IL}  |A_{II} ]   [u_I]  [G]
>
> and G should correspond to the Neumann problem. I might be thinking wrong
> (or maybe I completely misunderstood the idea), if so please correct me.
> But I think that the problem is that I am not explicitly telling PCBDDC
> which dofs are interface dofs.
>
> Regards,
> Abdullah Ali Sivas
>
> On Tue, 23 Oct 2018 at 23:16, Jed Brown  wrote:
>
>> Did you assemble "Neumann" problems that are compatible with your
>> definition of interior/interface degrees of freedom?
>>
>> Abdullah Ali Sivas  writes:
>>
>> > Dear all,
>> >
>> > I have a series of linear systems coming from a PDE for which BDDC is an
>> > optimal preconditioner. These linear systems are assembled and I read
>> them
>> > from a file, then convert into MATIS as required (as in
>> >
>> https://www.mcs.anl.gov/petsc/petsc-current/src/ksp/ksp/examples/tutorials/ex72.c.html
>> > ). I expect each of the systems converge to the solution in almost same
>> > number of iterations but I don't observe it. I think it is because I do
>> not
>> > provide enough information to the preconditioner. I can get a list of
>> inner
>> > dofs and interface dofs. However, I do not know how to use them. Has
>> anyone
>> > have any insights about it or done something similar?
>> >
>> > Best wishes,
>> > Abdullah Ali Sivas
>>
>

-- 
Stefano

Re: [petsc-users] [SLEPc] Number of iterations changes with MPI processes in Lanczos

[Jose E. Roman]

This is very strange. Make sure you call EPSSetFromOptions in the code. Do 
iteration counts change also for two different runs with the same number of 
processes?
Maybe Lanczos with default options is too sensitive (by default it does not 
reorthogonalize). Suggest using Krylov-Schur or Lanczos with full 
reorthogonalization (EPSLanczosSetReorthog).
Also, send the output of -eps_view to see if there is anything abnormal.

Jose


> El 24 oct 2018, a las 9:09, Ale Foggia  escribió:
> 
> I've tried the option that you gave me but I still get different number of 
> iterations when changing the number of MPI processes: I did 960 procs and 
> 1024 procs and I got 435 and 176 iterations, respectively.
> 
> El mar., 23 oct. 2018 a las 16:48, Jose E. Roman () 
> escribió:
> 
> 
> > El 23 oct 2018, a las 15:46, Ale Foggia  escribió:
> > 
> > 
> > 
> > El mar., 23 oct. 2018 a las 15:33, Jose E. Roman () 
> > escribió:
> > 
> > 
> > > El 23 oct 2018, a las 15:17, Ale Foggia  escribió:
> > > 
> > > Hello Jose, thanks for your answer.
> > > 
> > > El mar., 23 oct. 2018 a las 12:59, Jose E. Roman () 
> > > escribió:
> > > There is an undocumented option:
> > > 
> > >   -bv_reproducible_random
> > > 
> > > It will force the initial vector of the Krylov subspace to be the same 
> > > irrespective of the number of MPI processes. This should be used for 
> > > scaling analyses as the one you are trying to do.
> > > 
> > > What about when I'm not doing the scaling? Now I would like to ask for 
> > > computing time for bigger size problems, should I also use this option in 
> > > that case? Because, what happens if I have a "bad" configuration? 
> > > Meaning, I ask for some time, enough if I take into account the "correct" 
> > > scaling, but when I run it takes double the time/iterations, like it 
> > > happened before when changing from 960 to 1024 processes?
> > 
> > When you increase the matrix size the spectrum of the matrix changes and 
> > probably also the convergence, so the computation time is not easy to 
> > predict in advance.
> > 
> > Okey, I'll keep that in mine. I thought that, even if the spectrum changes, 
> > if I had a behaviour/tendency for 6 or 7 smaller cases I could predict more 
> > or less the time. It was working this way until I found this "iterations 
> > problem" which doubled the time of execution for the same size problem. To 
> > be completely sure, do you suggest me or not to use this run-time option 
> > when going in production? Can you elaborate a bit in the effect this 
> > option? Is the (huge) difference I got in the number of iterations 
> > something expected?
> 
> Ideally if you have a rough approximation of the eigenvector, you set it as 
> the initial vector with EPSSetInitialSpace(). Otherwise, SLEPc generates a 
> random initial vector, that is start the search blindly. The difference 
> between using one random vector or another may be large, depending on the 
> problem. Krylov-Schur is usually less sensitive to the initial vector. 
> 
> Jose
> 
> > 
> > 
> > > 
> > > An additional comment is that we strongly recommend to use the default 
> > > solver (Krylov-Schur), which will do Lanczos with implicit restart. It is 
> > > generally faster and more stable.
> > > 
> > > I will be doing Dynamical Lanczos, that means that I'll need the "matrix 
> > > whose rows are the eigenvectors of the tridiagonal matrix" (so, according 
> > > to the Lanczos Technical Report notation, I need the "matrix whose rows 
> > > are the eigenvectors of T_m", which should be the same as the vectors 
> > > y_i). I checked the Technical Report for Krylov-Schur also and I think I 
> > > can get the same information also from that solver, but I'm not sure. Can 
> > > you confirm this please? 
> > > Also, as the vectors I want are given by V_m^(-1)*x_i=y_i (following the 
> > > notation on the Report), my idea to get them was to retrieve the 
> > > invariant subspace V_m (with EPSGetInvariantSubspace), invert it, and 
> > > then multiply it with the eigenvectors that I get with EPSGetEigenvector. 
> > > Is there another easier (or with less computations) way to get this?
> > 
> > In Krylov-Schur the tridiagonal matrix T_m becomes 
> > arrowhead-plus-tridiagonal. Apart from this, it should be equivalent. The 
> > relevant information can be obtained with EPSGetBV() and EPSGetDS(). But 
> > this is a "developer level" interface. We could help you get this running. 
> > Send a small problem matrix to slepc-maint together with a more detailed 
> > description of what you need to compute.
> > 
> > Thanks! When I get to that part I'll write to slepc-maint for help.
> > 
> > 
> > Jose
> > 
> > > 
> > > 
> > > Jose
> > > 
> > > 
> > > 
> > > > El 23 oct 2018, a las 12:13, Ale Foggia  escribió:
> > > > 
> > > > Hello, 
> > > > 
> > > > I'm currently using Lanczos solver (EPSLANCZOS) to get the smallest 
> > > > real eigenvalue (EPS_SMALLEST_REAL) of a Hermitian problem (EPS_HEP). 
> > > > Those are the only options 

Re: [petsc-users] [SLEPc] Number of iterations changes with MPI processes in Lanczos

[Ale Foggia]

I've tried the option that you gave me but I still get different number of
iterations when changing the number of MPI processes: I did 960 procs and
1024 procs and I got 435 and 176 iterations, respectively.

El mar., 23 oct. 2018 a las 16:48, Jose E. Roman ()
escribió:

>
>
> > El 23 oct 2018, a las 15:46, Ale Foggia  escribió:
> >
> >
> >
> > El mar., 23 oct. 2018 a las 15:33, Jose E. Roman ()
> escribió:
> >
> >
> > > El 23 oct 2018, a las 15:17, Ale Foggia  escribió:
> > >
> > > Hello Jose, thanks for your answer.
> > >
> > > El mar., 23 oct. 2018 a las 12:59, Jose E. Roman ()
> escribió:
> > > There is an undocumented option:
> > >
> > >   -bv_reproducible_random
> > >
> > > It will force the initial vector of the Krylov subspace to be the same
> irrespective of the number of MPI processes. This should be used for
> scaling analyses as the one you are trying to do.
> > >
> > > What about when I'm not doing the scaling? Now I would like to ask for
> computing time for bigger size problems, should I also use this option in
> that case? Because, what happens if I have a "bad" configuration? Meaning,
> I ask for some time, enough if I take into account the "correct" scaling,
> but when I run it takes double the time/iterations, like it happened before
> when changing from 960 to 1024 processes?
> >
> > When you increase the matrix size the spectrum of the matrix changes and
> probably also the convergence, so the computation time is not easy to
> predict in advance.
> >
> > Okey, I'll keep that in mine. I thought that, even if the spectrum
> changes, if I had a behaviour/tendency for 6 or 7 smaller cases I could
> predict more or less the time. It was working this way until I found this
> "iterations problem" which doubled the time of execution for the same size
> problem. To be completely sure, do you suggest me or not to use this
> run-time option when going in production? Can you elaborate a bit in the
> effect this option? Is the (huge) difference I got in the number of
> iterations something expected?
>
> Ideally if you have a rough approximation of the eigenvector, you set it
> as the initial vector with EPSSetInitialSpace(). Otherwise, SLEPc generates
> a random initial vector, that is start the search blindly. The difference
> between using one random vector or another may be large, depending on the
> problem. Krylov-Schur is usually less sensitive to the initial vector.
>
> Jose
>
> >
> >
> > >
> > > An additional comment is that we strongly recommend to use the default
> solver (Krylov-Schur), which will do Lanczos with implicit restart. It is
> generally faster and more stable.
> > >
> > > I will be doing Dynamical Lanczos, that means that I'll need the
> "matrix whose rows are the eigenvectors of the tridiagonal matrix" (so,
> according to the Lanczos Technical Report notation, I need the "matrix
> whose rows are the eigenvectors of T_m", which should be the same as the
> vectors y_i). I checked the Technical Report for Krylov-Schur also and I
> think I can get the same information also from that solver, but I'm not
> sure. Can you confirm this please?
> > > Also, as the vectors I want are given by V_m^(-1)*x_i=y_i (following
> the notation on the Report), my idea to get them was to retrieve the
> invariant subspace V_m (with EPSGetInvariantSubspace), invert it, and then
> multiply it with the eigenvectors that I get with EPSGetEigenvector. Is
> there another easier (or with less computations) way to get this?
> >
> > In Krylov-Schur the tridiagonal matrix T_m becomes
> arrowhead-plus-tridiagonal. Apart from this, it should be equivalent. The
> relevant information can be obtained with EPSGetBV() and EPSGetDS(). But
> this is a "developer level" interface. We could help you get this running.
> Send a small problem matrix to slepc-maint together with a more detailed
> description of what you need to compute.
> >
> > Thanks! When I get to that part I'll write to slepc-maint for help.
> >
> >
> > Jose
> >
> > >
> > >
> > > Jose
> > >
> > >
> > >
> > > > El 23 oct 2018, a las 12:13, Ale Foggia 
> escribió:
> > > >
> > > > Hello,
> > > >
> > > > I'm currently using Lanczos solver (EPSLANCZOS) to get the smallest
> real eigenvalue (EPS_SMALLEST_REAL) of a Hermitian problem (EPS_HEP). Those
> are the only options I set for the solver. My aim is to be able to
> predict/estimate the time-to-solution. To do so, I was doing a scaling of
> the code for different sizes of matrices and for different number of MPI
> processes. As I was not observing a good scaling I checked the number of
> iterations of the solver (given by EPSGetIterationNumber). I've encounter
> that for the **same size** of matrix (that meaning, the same problem), when
> I change the number of MPI processes, the amount of iterations changes, and
> the behaviour is not monotonic. This are the numbers I've got:
> > > >
> > > > # procs   # iters
> > > > 960  157
> > > > 992  189
> > > > 1024338
> > > > 1056