Re: [petsc-users] DMPLEX project function error

[Matthew Knepley]

On Thu, Oct 11, 2018 at 4:39 PM Ellen M. Price 
wrote:

> I was working with a DMPLEX and FEM, following SNES example 12. I get
> the following error when I call DMProjectFunction, but I don't know what
> it means. Can anyone explain where I might have gone wrong, or at least
> what this error is telling me? I think the point closure size is
> correct, since my mesh is 3d simplex,


Yes, if you have 3D SIMPLEX mesh and are using P1 elements, then you would
have
4 dofs in the closure of a cell. The dual space dimension is the number of
dual space
basis vectors assigned to points in the closure. Since it is 1, it looks
like you have a P0
dual space. I assume you changed something in ex12?

  Thanks,

 Matt


> but what is the dual space
> dimension, and where might I have set it incorrectly?
>
> [0]PETSC ERROR: Nonconforming object sizes
> [0]PETSC ERROR: The section point closure size 4 != dual space dimension 1
> [0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html
> for trouble shooting.
> [0]PETSC ERROR: Petsc Release Version 3.9.2, May, 20, 2018
> ...
> [0]PETSC ERROR: Configure options
> --prefix=/home/eprice/software/petsc-opt --with-hdf5=1
> --with-hdf5-dir=/home/eprice/software/hdf5-parallel --with-mpe=1
> --with-mpe-dir=/home/eprice/software/mpe --with-debugging=0
> LDFLAGS="-pthread -lz" COPTFLAGS="-O3 -march=native -mtune=native"
> CXXOPTFLAGS="-O3 -march=native -mtune=native" FOPTFLAGS="-O3
> -march=native -mtune=native" --with-mpi=1
> --with-mpi-dir=/home/eprice/software/mpich --with-mumps=1
> --with-mumps-dir=/home/eprice/software/mumps --with-parmetis=1
> --with-parmetis-dir=/home/eprice/software/parmetis --with-metis=1
> --with-metis-dir=/home/eprice/software/parmetis --with-ptscotch=1
> --with-ptscotch-dir=/home/eprice/software/scotch --with-scalapack=1
> --with-scalapack-dir=/home/eprice/software/scalapack
> [0]PETSC ERROR: #1 DMProjectLocal_Generic_Plex() line 347 in
> /h/sabriel0/src/petsc-3.9.2/src/dm/impls/plex/plexproject.c
> [0]PETSC ERROR: #2 DMProjectFunctionLocal_Plex() line 428 in
> /h/sabriel0/src/petsc-3.9.2/src/dm/impls/plex/plexproject.c
> [0]PETSC ERROR: #3 DMProjectFunctionLocal() line 6265 in
> /h/sabriel0/src/petsc-3.9.2/src/dm/interface/dm.c
> [0]PETSC ERROR: #4 DMProjectFunction() line 6250 in
> /h/sabriel0/src/petsc-3.9.2/src/dm/interface/dm.c
> ...
>
> (I know this is an optimized PETSc build, but I get the same error from
> my debug build, it's just much slower.)
>
> Ellen
>


-- 
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] Failure of MUMPS

[Matthew Knepley]

On Thu, Oct 11, 2018 at 3:26 PM Michael Wick 
wrote:

> Thanks for all the suggestions!
>
> Increasing the value of icntl_14 in MUMPS helps a lot for my case.
>
> Do you have any suggestions for higher-order methods in saddle-point
> problems?
>

It depends on what the Schur complement looks like.

   Matt


> Mike
>
> Dave May  于2018年10月11日周四 上午1:50写道：
>
>>
>>
>> On Sat, 6 Oct 2018 at 12:42, Matthew Knepley  wrote:
>>
>>> On Fri, Oct 5, 2018 at 9:08 PM Mike Wick 
>>> wrote:
>>>
 Hello PETSc team:

 I am trying to solve a PDE problem with high-order finite elements. The
 matrix is getting denser and my experience is that MUMPS just outperforms
 iterative solvers.

>>>
>>> If the problem is elliptic, there is a lot of evidence that the P1
>>> preconditioner is descent for the system. Some people
>>> just project the system to P1, invert that with multigrid, and use that
>>> as the PC for Krylov. It should be worth trying.
>>>
>>
>> Matt means project to P1 directly from your high order function space in
>> one step. It is definitely worth trying.
>> For those interested, this approach is first described and discussed (to
>> my knowledge) in this paper:
>>
>> Persson, Per-Olof, and Jaime Peraire. "An efficient low memory implicit
>> DG algorithm for time dependent problems." *44th AIAA Aerospace Sciences
>> Meeting and Exhibit*. 2006.
>>
>>
>>> Moreover, as Jed will tell you, forming matrices for higher order is
>>> counterproductive. You should apply those matrix-free.
>>>
>>
>> I definitely agree with that.
>>
>> Cheers,
>>   Dave
>>
>>
>>
>>>
>>>   Thanks,
>>>
>>>  Matt
>>>
>>>
 For certain problems, MUMPS just fail in the middle for no clear
 reason. I just wander if there is any suggestion to improve the robustness
 of MUMPS? Or in general, any suggestion for interative solver with very
 high-order finite elements?

 Thanks!

 Mike

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

-- 
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] DMPLEX project function error

[Ellen M. Price]

I was working with a DMPLEX and FEM, following SNES example 12. I get
the following error when I call DMProjectFunction, but I don't know what
it means. Can anyone explain where I might have gone wrong, or at least
what this error is telling me? I think the point closure size is
correct, since my mesh is 3d simplex, but what is the dual space
dimension, and where might I have set it incorrectly?

[0]PETSC ERROR: Nonconforming object sizes
[0]PETSC ERROR: The section point closure size 4 != dual space dimension 1
[0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html
for trouble shooting.
[0]PETSC ERROR: Petsc Release Version 3.9.2, May, 20, 2018
...
[0]PETSC ERROR: Configure options
--prefix=/home/eprice/software/petsc-opt --with-hdf5=1
--with-hdf5-dir=/home/eprice/software/hdf5-parallel --with-mpe=1
--with-mpe-dir=/home/eprice/software/mpe --with-debugging=0
LDFLAGS="-pthread -lz" COPTFLAGS="-O3 -march=native -mtune=native"
CXXOPTFLAGS="-O3 -march=native -mtune=native" FOPTFLAGS="-O3
-march=native -mtune=native" --with-mpi=1
--with-mpi-dir=/home/eprice/software/mpich --with-mumps=1
--with-mumps-dir=/home/eprice/software/mumps --with-parmetis=1
--with-parmetis-dir=/home/eprice/software/parmetis --with-metis=1
--with-metis-dir=/home/eprice/software/parmetis --with-ptscotch=1
--with-ptscotch-dir=/home/eprice/software/scotch --with-scalapack=1
--with-scalapack-dir=/home/eprice/software/scalapack
[0]PETSC ERROR: #1 DMProjectLocal_Generic_Plex() line 347 in
/h/sabriel0/src/petsc-3.9.2/src/dm/impls/plex/plexproject.c
[0]PETSC ERROR: #2 DMProjectFunctionLocal_Plex() line 428 in
/h/sabriel0/src/petsc-3.9.2/src/dm/impls/plex/plexproject.c
[0]PETSC ERROR: #3 DMProjectFunctionLocal() line 6265 in
/h/sabriel0/src/petsc-3.9.2/src/dm/interface/dm.c
[0]PETSC ERROR: #4 DMProjectFunction() line 6250 in
/h/sabriel0/src/petsc-3.9.2/src/dm/interface/dm.c
...

(I know this is an optimized PETSc build, but I get the same error from
my debug build, it's just much slower.)

Ellen

Re: [petsc-users] Failure of MUMPS

[Michael Wick]

Thanks for all the suggestions!

Increasing the value of icntl_14 in MUMPS helps a lot for my case.

Do you have any suggestions for higher-order methods in saddle-point
problems?

Mike

Dave May  于2018年10月11日周四 上午1:50写道：

>
>
> On Sat, 6 Oct 2018 at 12:42, Matthew Knepley  wrote:
>
>> On Fri, Oct 5, 2018 at 9:08 PM Mike Wick 
>> wrote:
>>
>>> Hello PETSc team:
>>>
>>> I am trying to solve a PDE problem with high-order finite elements. The
>>> matrix is getting denser and my experience is that MUMPS just outperforms
>>> iterative solvers.
>>>
>>
>> If the problem is elliptic, there is a lot of evidence that the P1
>> preconditioner is descent for the system. Some people
>> just project the system to P1, invert that with multigrid, and use that
>> as the PC for Krylov. It should be worth trying.
>>
>
> Matt means project to P1 directly from your high order function space in
> one step. It is definitely worth trying.
> For those interested, this approach is first described and discussed (to
> my knowledge) in this paper:
>
> Persson, Per-Olof, and Jaime Peraire. "An efficient low memory implicit DG
> algorithm for time dependent problems." *44th AIAA Aerospace Sciences
> Meeting and Exhibit*. 2006.
>
>
>> Moreover, as Jed will tell you, forming matrices for higher order is
>> counterproductive. You should apply those matrix-free.
>>
>
> I definitely agree with that.
>
> Cheers,
>   Dave
>
>
>
>>
>>   Thanks,
>>
>>  Matt
>>
>>
>>> For certain problems, MUMPS just fail in the middle for no clear reason.
>>> I just wander if there is any suggestion to improve the robustness of
>>> MUMPS? Or in general, any suggestion for interative solver with very
>>> high-order finite elements?
>>>
>>> Thanks!
>>>
>>> Mike
>>>
>>
>>
>> --
>> 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] STFILTER in slepc

[Jose E. Roman]

Run your code with -eps_view_mat0 binary:amatrix.bin
This will save the A matrix of the EPS solver in file "amatrix.bin" at the end 
of EPSSolve().
I cannot say which is the best method, it depends on the distribution of the 
spectrum. Let's see how it goes.

Jose


> El 11 oct 2018, a las 20:13, Moritz Cygorek  escribió:
> 
> Thank you for the response. This example is exactly the kind of thing I 
> thought about.
> But, as you say, convergence is indeed not better. 
> 
> 
> What's the best format to send the matrix to slepc-maint? 
> My small test case is a sparse matrix of dimension 30720. I currently read it 
> from three files (containing IA, JA, and the values in CSR format). Because 
> everything is contained in a larger framework, it's not very easy to extract 
> the reader.  Is there a generic way to read in sparse matrices from files for 
> PETSc?
> 
> 
> 
> Maybe you can already give me a hint about the best method if I tell you the 
> properties of the spectrum:
> The problem is the calculation of electronic states in a quantum dot.  There 
> is a band gap from around 0 to 2 and I am interested in the first 20 or 40 
> eigenvectors above and below the gap.
> The first few states are separated by about 0.1 but for higher states the 
> energies come closer and closer together.
> Additionally, there are a number of states with energies around 1000, which 
> are artificial and originate from the way we treat the boundary conditions. 
> Also, I know that all states come in pairs with the same energy (Kramers 
> degeneracy).
> 
> For the small test case, the conduction band states, i.e. eigenvectors with 
> energies close to 2, converge very fast (about 5 min on a laptop computer). 
> However, the states with energies around 0 converge much more slowly and 
> that's one of my major problems. 
> For those states harmonic extraction seems to be better suited, but I have 
> the impression that it is not extremely stable. For example, applied to the 
> states close to 2 I see that some states are skipped, which can be seen by 
> the fact that the degeneracies are sometimes wrong. 
> Also, with harmonic extraction, the program sometimes stops claiming the 
> number of requested eigenpairs are reached, but the calculated relative error 
> of most states is way larger than the tolerance.
> 
> Maybe you know from experience which method is better suited to tackle these 
> kinds of problems?
> 
> 
> Eventually, I intend to do calculates with dimensions of ~ 10 million 
> distributed of a few 100 CPUs.
> 
> Regards,
> Moritz
> 
> 
> 
> 
> From: Jose E. Roman 
> Sent: Thursday, October 11, 2018 5:55 AM
> To: Matthew Knepley; Moritz Cygorek
> Cc: PETSc; Carmen Campos
> Subject: Re: [petsc-users] STFILTER in slepc
>  
> The filter technique must be used with an interval in EPS, e.g.
>   -eps_interval 2.,2.7
> (This interval will be passed to STFILTER, so no need to specify 
> -st_filter_interval).
> Therefore, it does not require -eps_target_magnitude or similar.
> 
> Regarding the simpler (A-tau*I)^2, we call it "spectrum folding" (see section 
> 3.4.6 of SLEPc's manual) and it is implemented in ex24.c
> http://slepc.upv.es/documentation/current/src/eps/examples/tutorials/ex24.c.html
> 
> My experience is that spectrum folding will not give good convergence except 
> for easy problems. For more difficult problems with clustered eigenvalues you 
> need a polynomial filter. The technique of polynomial filter trades 
> iterations of the Krylov subspace for matrix-vector products required by a 
> high-degree polynomial.
> 
> If you want, send us the matrix to slepc-maint and we will have a try.
> 
> Jose
>  
> 
> > El 11 oct 2018, a las 2:43, Matthew Knepley  escribió:
> > 
> > On Wed, Oct 10, 2018 at 8:41 PM Moritz Cygorek  wrote:
> > Thank you very much. Apparently, I've misunderstood what the filter 
> > actually does. I thought about the much simpler process, where you 
> > diagonalize
> > 
> > 
> > 
> > -(A- tau*I)^2 +offset*I
> > 
> > 
> > 
> > where tau is my target an offset is large enough so that the global maximum 
> > is reached for eigenvalues around tau.
> > 
> > 
> > Is this different from -eps_target_magnitude?
> > 
> >   Thanks,
> > 
> > Matt
> >  
> > Then you look for the largest eigenvalue of the modified problem and either 
> > calculate the Ritz value of the original matrix or calculate back from the 
> > eigenvalues of the modified problem.
> > 
> > 
> > 
> > Now, it looks to me like -st_type filter activates something like the 
> > package FILTLAN.
> > 
> > 
> > 
> > I guess I can define a MatShell to do the thing I intended in the first 
> > place.
> > 
> > But, I guess, this is a common thing, so I am wondering whether it is 
> > already implemented somewhere and I just didn't find it in the 
> > documentation.  Can you say something about this?
> > 
> > 
> > 
> > Regards,
> > 
> > Moritz
> > 
> > 
> > 
> > 
> > 
> > From: Jose E. Roman 
> > Sent: Wednesday, October 10, 2018 3:48 PM
> 

Re: [petsc-users] STFILTER in slepc

[Moritz Cygorek]

Thank you for the response. This example is exactly the kind of thing I thought 
about.

But, as you say, convergence is indeed not better.



What's the best format to send the matrix to slepc-maint?

My small test case is a sparse matrix of dimension 30720. I currently read it 
from three files (containing IA, JA, and the values in CSR format). Because 
everything is contained in a larger framework, it's not very easy to extract 
the reader.  Is there a generic way to read in sparse matrices from files for 
PETSc?




Maybe you can already give me a hint about the best method if I tell you the 
properties of the spectrum:

The problem is the calculation of electronic states in a quantum dot.  There is 
a band gap from around 0 to 2 and I am interested in the first 20 or 40 
eigenvectors above and below the gap.

The first few states are separated by about 0.1 but for higher states the 
energies come closer and closer together.

Additionally, there are a number of states with energies around 1000, which are 
artificial and originate from the way we treat the boundary conditions. Also, I 
know that all states come in pairs with the same energy (Kramers degeneracy).


For the small test case, the conduction band states, i.e. eigenvectors with 
energies close to 2, converge very fast (about 5 min on a laptop computer). 
However, the states with energies around 0 converge much more slowly and that's 
one of my major problems.

For those states harmonic extraction seems to be better suited, but I have the 
impression that it is not extremely stable. For example, applied to the states 
close to 2 I see that some states are skipped, which can be seen by the fact 
that the degeneracies are sometimes wrong.

Also, with harmonic extraction, the program sometimes stops claiming the number 
of requested eigenpairs are reached, but the calculated relative error of most 
states is way larger than the tolerance.


Maybe you know from experience which method is better suited to tackle these 
kinds of problems?



Eventually, I intend to do calculates with dimensions of ~ 10 million 
distributed of a few 100 CPUs.


Regards,
Moritz





From: Jose E. Roman 
Sent: Thursday, October 11, 2018 5:55 AM
To: Matthew Knepley; Moritz Cygorek
Cc: PETSc; Carmen Campos
Subject: Re: [petsc-users] STFILTER in slepc

The filter technique must be used with an interval in EPS, e.g.
  -eps_interval 2.,2.7
(This interval will be passed to STFILTER, so no need to specify 
-st_filter_interval).
Therefore, it does not require -eps_target_magnitude or similar.

Regarding the simpler (A-tau*I)^2, we call it "spectrum folding" (see section 
3.4.6 of SLEPc's manual) and it is implemented in ex24.c
http://slepc.upv.es/documentation/current/src/eps/examples/tutorials/ex24.c.html

My experience is that spectrum folding will not give good convergence except 
for easy problems. For more difficult problems with clustered eigenvalues you 
need a polynomial filter. The technique of polynomial filter trades iterations 
of the Krylov subspace for matrix-vector products required by a high-degree 
polynomial.

If you want, send us the matrix to slepc-maint and we will have a try.

Jose


> El 11 oct 2018, a las 2:43, Matthew Knepley  escribió:
>
> On Wed, Oct 10, 2018 at 8:41 PM Moritz Cygorek  wrote:
> Thank you very much. Apparently, I've misunderstood what the filter actually 
> does. I thought about the much simpler process, where you diagonalize
>
>
>
> -(A- tau*I)^2 +offset*I
>
>
>
> where tau is my target an offset is large enough so that the global maximum 
> is reached for eigenvalues around tau.
>
>
> Is this different from -eps_target_magnitude?
>
>   Thanks,
>
> Matt
>
> Then you look for the largest eigenvalue of the modified problem and either 
> calculate the Ritz value of the original matrix or calculate back from the 
> eigenvalues of the modified problem.
>
>
>
> Now, it looks to me like -st_type filter activates something like the package 
> FILTLAN.
>
>
>
> I guess I can define a MatShell to do the thing I intended in the first place.
>
> But, I guess, this is a common thing, so I am wondering whether it is already 
> implemented somewhere and I just didn't find it in the documentation.  Can 
> you say something about this?
>
>
>
> Regards,
>
> Moritz
>
>
>
>
>
> From: Jose E. Roman 
> Sent: Wednesday, October 10, 2018 3:48 PM
> To: Moritz Cygorek
> Cc: petsc-users@mcs.anl.gov
> Subject: Re: [petsc-users] STFILTER in slepc
>
> This type of method requires a very high degree polynomial; suggest using 
> degree=100 at least (this is the default value), but larger values may be 
> necessary. Also, for this particular filter the "range" must be approximately 
> equal to the numerical range; if you have no clue where your first and last 
> eigenvalues are, you may use EPSSolve() calls with EPS_LARGEST_REAL and 
> EPS_SMALLEST_REAL.
>
> Jose
>
> > El 10 oct 2018, a las 21:10, Moritz Cygorek  escribió:
> >
> > 

Re: [petsc-users] Problems about Matrix-Free and DM object

[Matthew Knepley]

On Thu, Oct 11, 2018 at 11:38 AM Yingjie Wu  wrote:

> Dear Petsc developer:
> Hi,
> Thank you very much for your previous reply. I feel that I have learned a
> lot about the matrix-free method.
> I was studying the example /src/snes/example/tutorials/ex28.c recently.
> One of its tests is as follows:
>
> -u_da_grid_x 20 -snes_converged_reason -snes_monitor_short
> -ksp_monitor_short -problem_type 2 -snes_mf_operator -pack_dm_mat_type
> {{aij nest}} -pc_type fieldsplit -pc_fieldsplit_dm_splits
> -pc_fieldsplit_type additive -fieldsplit_u_ksp_type gmres
> -fieldsplit_k_pc_type jacobi
>
> The example is calculated by the Matrix-free method, and the precondition
> matrix is divided by FieldSplit. There are several questions in the
> procedure that make me feel uncertain.
> 1. The program uses the Matrixfree method by the command
> "-snes_mf_operator", but the form of Jacobian matrix is not defined as
> MATMFFD in the program. Does the program automatically set the Jacobian
> matrix to the form of MATMFFD after opening the command?
>

If you give that option, PETSc will automatically create a MFFD Mat object
and stick your nonlinear residual evaluation in it.


> Both Jacobian matrix and precondition matrix are defined as B in the
> program, and matrix elements are provided by Jacobian evaluation program.
> But the Jacobian matrix in the Matrixfree method is in the form of MATMFFD.
> You mentioned earlier that adding elements to the matrix with type MATMFFD
> can cause program errors, but why does the program work in this example?
> The procedures are as follows:
>
> ierr = SNESSetFunction(snes,F,FormFunction_All,user);CHKERRQ(ierr);
> ierr = SNESSetJacobian(snes,B,B,FormJacobian_All,user);CHKERRQ(ierr);
>
>
Because the FoirmJacobian() function in ex28 only ever puts values into the
preconditioning matrix. This is generally the right way to proceed.


> 2.Can we use DMCompositeScatter to scatter ghost value into local vectors?
> I didn't seem to mention the documentation. The program only uses the
> following functions, but can use ghost value.
> ierr = DMCompositeScatter(user->pack,X,Uloc,Kloc);CHKERRQ(ierr);
> ierr = DMDAVecGetArray(dau,Uloc,);CHKERRQ(ierr);
> ierr = DMDAVecGetArray(dak,Kloc,);CHKERRQ(ierr);
>

DMCompositeScatter splits a vector up into the smaller vectors for each
system. DMComposite is very specialized, and we do not
recommend using it except in special cases.

DMGlobalToLocal() and DMLocalToGlobal() communicate ghost values.

  Thanks,

 Matt


> Greatly appreciated your reply and attention.
> 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 about Matrix-Free and DM object

[Yingjie Wu]

Dear Petsc developer:
Hi,
Thank you very much for your previous reply. I feel that I have learned a
lot about the matrix-free method.
I was studying the example /src/snes/example/tutorials/ex28.c recently. One
of its tests is as follows:

-u_da_grid_x 20 -snes_converged_reason -snes_monitor_short
-ksp_monitor_short -problem_type 2 -snes_mf_operator -pack_dm_mat_type
{{aij nest}} -pc_type fieldsplit -pc_fieldsplit_dm_splits
-pc_fieldsplit_type additive -fieldsplit_u_ksp_type gmres
-fieldsplit_k_pc_type jacobi

The example is calculated by the Matrix-free method, and the precondition
matrix is divided by FieldSplit. There are several questions in the
procedure that make me feel uncertain.
1. The program uses the Matrixfree method by the command
"-snes_mf_operator", but the form of Jacobian matrix is not defined as
MATMFFD in the program. Does the program automatically set the Jacobian
matrix to the form of MATMFFD after opening the command?

Both Jacobian matrix and precondition matrix are defined as B in the
program, and matrix elements are provided by Jacobian evaluation program.
But the Jacobian matrix in the Matrixfree method is in the form of MATMFFD.
You mentioned earlier that adding elements to the matrix with type MATMFFD
can cause program errors, but why does the program work in this example?
The procedures are as follows:

ierr = SNESSetFunction(snes,F,FormFunction_All,user);CHKERRQ(ierr);
ierr = SNESSetJacobian(snes,B,B,FormJacobian_All,user);CHKERRQ(ierr);


2.Can we use DMCompositeScatter to scatter ghost value into local vectors?
I didn't seem to mention the documentation. The program only uses the
following functions, but can use ghost value.
ierr = DMCompositeScatter(user->pack,X,Uloc,Kloc);CHKERRQ(ierr);
ierr = DMDAVecGetArray(dau,Uloc,);CHKERRQ(ierr);
ierr = DMDAVecGetArray(dak,Kloc,);CHKERRQ(ierr);

Greatly appreciated your reply and attention.
Thanks,
Yingjie

Re: [petsc-users] STFILTER in slepc

[Jose E. Roman]

The filter technique must be used with an interval in EPS, e.g.
  -eps_interval 2.,2.7
(This interval will be passed to STFILTER, so no need to specify 
-st_filter_interval).
Therefore, it does not require -eps_target_magnitude or similar.

Regarding the simpler (A-tau*I)^2, we call it "spectrum folding" (see section 
3.4.6 of SLEPc's manual) and it is implemented in ex24.c
http://slepc.upv.es/documentation/current/src/eps/examples/tutorials/ex24.c.html

My experience is that spectrum folding will not give good convergence except 
for easy problems. For more difficult problems with clustered eigenvalues you 
need a polynomial filter. The technique of polynomial filter trades iterations 
of the Krylov subspace for matrix-vector products required by a high-degree 
polynomial.

If you want, send us the matrix to slepc-maint and we will have a try.

Jose
 

> El 11 oct 2018, a las 2:43, Matthew Knepley  escribió:
> 
> On Wed, Oct 10, 2018 at 8:41 PM Moritz Cygorek  wrote:
> Thank you very much. Apparently, I've misunderstood what the filter actually 
> does. I thought about the much simpler process, where you diagonalize
> 
> 
> 
> -(A- tau*I)^2 +offset*I
> 
> 
> 
> where tau is my target an offset is large enough so that the global maximum 
> is reached for eigenvalues around tau.
> 
> 
> Is this different from -eps_target_magnitude?
> 
>   Thanks,
> 
> Matt
>  
> Then you look for the largest eigenvalue of the modified problem and either 
> calculate the Ritz value of the original matrix or calculate back from the 
> eigenvalues of the modified problem.
> 
> 
> 
> Now, it looks to me like -st_type filter activates something like the package 
> FILTLAN.
> 
> 
> 
> I guess I can define a MatShell to do the thing I intended in the first place.
> 
> But, I guess, this is a common thing, so I am wondering whether it is already 
> implemented somewhere and I just didn't find it in the documentation.  Can 
> you say something about this?
> 
> 
> 
> Regards,
> 
> Moritz
> 
> 
> 
> 
> 
> From: Jose E. Roman 
> Sent: Wednesday, October 10, 2018 3:48 PM
> To: Moritz Cygorek
> Cc: petsc-users@mcs.anl.gov
> Subject: Re: [petsc-users] STFILTER in slepc
>  
> This type of method requires a very high degree polynomial; suggest using 
> degree=100 at least (this is the default value), but larger values may be 
> necessary. Also, for this particular filter the "range" must be approximately 
> equal to the numerical range; if you have no clue where your first and last 
> eigenvalues are, you may use EPSSolve() calls with EPS_LARGEST_REAL and 
> EPS_SMALLEST_REAL.
> 
> Jose
> 
> > El 10 oct 2018, a las 21:10, Moritz Cygorek  escribió:
> > 
> > Thank you for the fast reply. 
> > 
> > I've tried running my program (using the defaul Krylov-Schur method for 
> > sparse MPI matrices) with the additional options:
> > 
> > -st_type filter -st_filter_degree 2 -st_filter_interval 2.,2.7  
> > -st_filter_range -2000,2000
> > 
> > and I get the following error message:
> > 
> > [0]PETSC ERROR: STFILTER cannot get the filter specified; please adjust 
> > your filter parameters (e.g. increasing the polynomial degree)
> > 
> > [0]PETSC ERROR: #1 FILTLAN_GetIntervals() line 451 in 
> > /home/applications/sources/libraries/slepc-3.9.2/src/sys/classes/st/impls/filter/filtlan.c
> > [0]PETSC ERROR: #2 STFilter_FILTLAN_setFilter() line 1016 in 
> > /home/applications/sources/libraries/slepc-3.9.2/src/sys/classes/st/impls/filter/filtlan.c
> > [0]PETSC ERROR: #3 STSetUp_Filter() line 42 in 
> > /home/applications/sources/libraries/slepc-3.9.2/src/sys/classes/st/impls/filter/filter.c
> > [0]PETSC ERROR: #4 STSetUp() line 271 in 
> > /home/applications/sources/libraries/slepc-3.9.2/src/sys/classes/st/interface/stsolve.c
> > [0]PETSC ERROR: #5 EPSSetUp() line 263 in 
> > /home/applications/sources/libraries/slepc-3.9.2/src/eps/interface/epssetup.c
> > [0]PETSC ERROR: #6 EPSSolve() line 135 in 
> > /home/applications/sources/libraries/slepc-3.9.2/src/eps/interface/epssolve.c
> > 
> > 
> > 
> > Do you have a clue what I've missed?
> > 
> > 
> > Moritz
> > 
> > 
> > From: Jose E. Roman 
> > Sent: Wednesday, October 10, 2018 2:30 PM
> > To: Moritz Cygorek
> > Cc: petsc-users@mcs.anl.gov
> > Subject: Re: [petsc-users] STFILTER in slepc
> >  
> > 
> > > El 10 oct 2018, a las 19:54, Moritz Cygorek  
> > > escribió:
> > > 
> > > Hi,
> > > 
> > > in the list of changes to SLEPc version 3.8, it is stated that there is a 
> > > preliminary implementation of polynomial filtering using STFILTER. 
> > > 
> > > Because I am struggling to obtain interior eigenvalues and harmonic 
> > > extraction seems not to be stable enough in my case, I wanted to give it 
> > > a try, but I could not find any documentation yet.
> > > 
> > > Does anybody have an example of how to use STFILTER or any documentation 
> > > about it?
> > > 
> > > Thanks in advance,
> > > Moritz
> > 
> > There are no examples. You just set the type to STFILTER and set some 
> > parameters such 

Re: [petsc-users] Failure of MUMPS

[Dave May]

On Sat, 6 Oct 2018 at 12:42, Matthew Knepley  wrote:

> On Fri, Oct 5, 2018 at 9:08 PM Mike Wick 
> wrote:
>
>> Hello PETSc team:
>>
>> I am trying to solve a PDE problem with high-order finite elements. The
>> matrix is getting denser and my experience is that MUMPS just outperforms
>> iterative solvers.
>>
>
> If the problem is elliptic, there is a lot of evidence that the P1
> preconditioner is descent for the system. Some people
> just project the system to P1, invert that with multigrid, and use that as
> the PC for Krylov. It should be worth trying.
>

Matt means project to P1 directly from your high order function space in
one step. It is definitely worth trying.
For those interested, this approach is first described and discussed (to my
knowledge) in this paper:

Persson, Per-Olof, and Jaime Peraire. "An efficient low memory implicit DG
algorithm for time dependent problems." *44th AIAA Aerospace Sciences
Meeting and Exhibit*. 2006.


> Moreover, as Jed will tell you, forming matrices for higher order is
> counterproductive. You should apply those matrix-free.
>

I definitely agree with that.

Cheers,
  Dave



>
>   Thanks,
>
>  Matt
>
>
>> For certain problems, MUMPS just fail in the middle for no clear reason.
>> I just wander if there is any suggestion to improve the robustness of
>> MUMPS? Or in general, any suggestion for interative solver with very
>> high-order finite elements?
>>
>> Thanks!
>>
>> Mike
>>
>
>
> --
> 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/
> 
>