Re: [petsc-users] Near nullspace lost in fieldsplit

[Jeremy Theler (External)]

>   The new code is in https://gitlab.com/petsc/petsc/-/merge_requests/7293 and 
> retains the null space on the submatrices for both MatZeroRows() and 
> MatZeroRowsAndColumns() regardless of changes to the nonzero structure of the 
> matrix.

Thank you very much Barry.
It works as expected.

--
jeremy

Re: [petsc-users] Near nullspace lost in fieldsplit

[Jeremy Theler (External)]

Hi Barry

>   The bug fix for 2 is availabel in 
> https://gitlab.com/petsc/petsc/-/merge_requests/7279

Note that our code goes through   MatZeroRowsColumns_MPIAIJ() not through 
MatZeroRows_MPIAIJ() so this fix does not change anything for the case I 
mentioned.

--
jeremy

Re: [petsc-users] Near nullspace lost in fieldsplit

[Jeremy Theler (External)]

> > Because of a combination of settings, our code passes through this line:
> >
> > https://gitlab.com/petsc/petsc/-/blob/main/src/ksp/pc/impls/fieldsplit/fieldsplit.c?ref_type=heads#L692
> >
> > i.e. the matrices associated with each of the sub-KSPs of a fieldsplit are 
> > destroyed and then re-created later.
> > The thing is that one of these destroyed matrices had a near nullspace 
> > attached, which is lost because the new matrix does > not have it anymore.
> >
> > Is this a bug or are we missing something?
>
> I just want to get a clear picture. You create a PCFIELDSPLIT, set it up, 
> then pull out the matrices and attach a nullspace before the solve.

We need to solve an SNES. We use dmplex so we have the jacobian allocated 
before starting the solve.
At setup time we

 1. define the PC of the KSP of the SNES to be fieldsplit
 2. define the fields with ISes
 3. call PCSetup() to create the sub-KSPs
 4. retrieve the matrix attached to the sub-KSP that needs the near nullspace 
and attach it to that matrix

> At a later time, you start another solve with this PC, and it has the 
> DIFFERENT_NONZERO_PATTERN flag, so it recreates these matrices and loses your 
> attached nullspace.

At a later time, in the jacobian evaluation we populate the global matrix (i.e. 
not the matrices attached to each sub-KSPs) and then we set dirichlet bcs with 
MatZeroRowsColumns() on that same global matrix.
For some reason, in serial the near nullspace is not lost but in parallel the 
call to MatZeroRowsColumns() does change the non-zero structure (even though 
the manual says it does not) and then the code goes through that line 692 in 
fieldsplit.c and the near nullspace is lost.

> First, does the matrix really change?

Well, the matrix during setup is not filled in, just allocated.
The thing is that if we set MAT_KEEP_NONZERO_PATTERN to true with 
MatSetOption() before setting the dirichlet BCs, then the near nullspace is not 
lost (because the code does not go through line 692 of fieldsplit.c).

So there are (at least) two issues:

 1. Code going through line 692 looses the near nullspace of the matrices 
attached to the sub-KSPs
 2. The call to MatZeroRowsColumns() changes then non-zero structure for MPIAIJ 
but not for SEQAIJ


> Second, I had the same problem. I added callbacks to DM which allow a 
> nullspace to be automatically attached if you extract a certain subfield. Are 
> you using a DM?


Yes. Can you give us an example?

Regards
--
jeremy

[petsc-users] Near nullspace lost in fieldsplit

[Jeremy Theler (External)]

Hi all

Because of a combination of settings, our code passes through this line:

https://gitlab.com/petsc/petsc/-/blob/main/src/ksp/pc/impls/fieldsplit/fieldsplit.c?ref_type=heads#L692

i.e. the matrices associated with each of the sub-KSPs of a fieldsplit are 
destroyed and then re-created later.
The thing is that one of these destroyed matrices had a near nullspace 
attached, which is lost because the new matrix does not have it anymore.

Is this a bug or are we missing something?
I'll need some time to come up with a MWE showing how the near nullspace is 
lost, but if its needed I can take a stab at it.

--
jeremy

[https://gitlab.com/uploads/-/system/project/avatar/13882401/PETSc_RBG-logo.png]
src/ksp/pc/impls/fieldsplit/fieldsplit.c · main · PETSc / petsc · 
GitLab
PETSc, pronounced PET-see (the S is silent), is a suite of data structures and 
routines for the scalable (parallel) solution of scientific applications 
modeled by partial differential equations.
gitlab.com

Re: [petsc-users] Help for MatNullSpaceCreateRigidBody

[Jeremy Theler (External)]

just in case it helps, here's one way I have to create the near nullspace:

https://github.com/seamplex/feenox/blob/main/src/pdes/mechanical/init.c#L468

--
jeremy

From: petsc-users  on behalf of Jordi Manyer 
Fuertes via petsc-users 
Sent: Tuesday, December 5, 2023 9:57 AM
To: Matthew Knepley ; bsm...@petsc.dev 
Cc: petsc-users@mcs.anl.gov 
Subject: Re: [petsc-users] Help for MatNullSpaceCreateRigidBody


[External Sender]

Thanks for the prompt response. Both answers look like what I'm doing.

After playing a bit more with solver, I managed to make it run in parallel with 
different boundary conditions (full dirichlet bcs, vs mixed newmann + 
dirichlet). This raises two questions:

- How relevant are boundary conditions (eliminating dirichlet rows/cols vs weak 
newmann bcs) to the solver? Should I modify something when changing boundary 
conditions?

- Also, the solver did well with the old bcs when run in a single processor 
(but not in parallel). This seems odd, since parallel and serial behavior 
should be consistent (or not?). Could it be fault of the PCGAMG? I believe the 
default local solver is ILU, shoud I be changing it to LU or something else for 
these kind of problems?

Thank you both again,

Jordi

Re: [petsc-users] Start logging with -info after PetscInitialize()

[Jeremy Theler (External) via petsc-users]

Thanks. That worked.

--
jeremy

From: Jacob Faibussowitsch 
Sent: Tuesday, June 13, 2023 11:25 AM
To: Jeremy Theler (External) 
Cc: petsc-users 
Subject: Re: [petsc-users] Start logging with -info after PetscInitialize()

[External Sender]

Call PetscInfoDestroy() first.

https://petsc.org/main/manualpages/Profiling/PetscInfoDestroy/

Best regards,

Jacob Faibussowitsch
(Jacob Fai - booss - oh - vitch)

> On Jun 13, 2023, at 10:16, Jeremy Theler (External) via petsc-users 
>  wrote:
>
> Hello all.
>
> I've asked this question to Satish personally last week at the conference, 
> but I'm stuck so any help would be appreciated.
> For some reason not worth explaining, I need to activate -info after 
> PetscInitialize() has been already called.
> I'm trying something like this:
>
> PetscOptionsSetValue(NULL, "-info", NULL);
> PetscInfoSetFromOptions(NULL);
>
> The second call fails with
>
> [0]PETSC ERROR: - Error Message 
> --
> [0]PETSC ERROR: Object is in wrong state
> [0]PETSC ERROR: PetscInfoSetClasses() cannot be called after 
> PetscInfoGetClass() or PetscInfoProcessClass()
> [0]PETSC ERROR: See https://petsc.org/release/faq/ for trouble shooting.
> [0]PETSC ERROR: Petsc Release Version 3.19.2, Jun 01, 2023
> [0]PETSC ERROR: reflexCLI on a double-int32-release named LIN54Z7SQ3 by 
> jtheler Tue Jun 13 11:08:29 2023
> [0]PETSC ERROR: Configure options --download-eigen --download-hdf5 
> --download-hypre --download-metis --download-mumps --download-parmetis 
> --download-pragmatic --download-scalapack --download-slepc 
> --with-64-bit-indices=no --with-debugging=no --with-precision=double 
> --with-scalar-type=real COPTFLAGS=-O3 CXXOPTFLAGS=-O3 FOPTFLAGS=-O3 
> --download-egads --download-opencascade --download-tetgen
> [0]PETSC ERROR: #1 PetscInfoSetClasses() at 
> /home/jtheler/libs/petsc-3.19.2/src/sys/info/verboseinfo.c:182
> [0]PETSC ERROR: #2 PetscInfoSetFromOptions() at 
> /home/jtheler/libs/petsc-3.19.2/src/sys/info/verboseinfo.c:407
>
> But if I ignore the non-zero return value and I allow my program to continue, 
> the required logging is enabled.
> I also tried using a local PetscOptions object but the result is the same.
>
> Any ideas to avoid that wrong state error?
>
> Thanks
> --
> jeremy

[petsc-users] Start logging with -info after PetscInitialize()

[Jeremy Theler (External) via petsc-users]

Hello all.

I've asked this question to Satish personally last week at the conference, but 
I'm stuck so any help would be appreciated.
For some reason not worth explaining, I need to activate -info after 
PetscInitialize() has been already called.
I'm trying something like this:

PetscOptionsSetValue(NULL, "-info", NULL);
PetscInfoSetFromOptions(NULL);

The second call fails with

[0]PETSC ERROR: - Error Message 
--
[0]PETSC ERROR: Object is in wrong state
[0]PETSC ERROR: PetscInfoSetClasses() cannot be called after 
PetscInfoGetClass() or PetscInfoProcessClass()
[0]PETSC ERROR: See https://petsc.org/release/faq/ for trouble shooting.
[0]PETSC ERROR: Petsc Release Version 3.19.2, Jun 01, 2023
[0]PETSC ERROR: reflexCLI on a double-int32-release named LIN54Z7SQ3 by jtheler 
Tue Jun 13 11:08:29 2023
[0]PETSC ERROR: Configure options --download-eigen --download-hdf5 
--download-hypre --download-metis --download-mumps --download-parmetis 
--download-pragmatic --download-scalapack --download-slepc 
--with-64-bit-indices=no --with-debugging=no --with-precision=double 
--with-scalar-type=real COPTFLAGS=-O3 CXXOPTFLAGS=-O3 FOPTFLAGS=-O3 
--download-egads --download-opencascade --download-tetgen
[0]PETSC ERROR: #1 PetscInfoSetClasses() at 
/home/jtheler/libs/petsc-3.19.2/src/sys/info/verboseinfo.c:182
[0]PETSC ERROR: #2 PetscInfoSetFromOptions() at 
/home/jtheler/libs/petsc-3.19.2/src/sys/info/verboseinfo.c:407

But if I ignore the non-zero return value and I allow my program to continue, 
the required logging is enabled.
I also tried using a local PetscOptions object but the result is the same.

Any ideas to avoid that wrong state error?

Thanks
--
jeremy

Re: [petsc-users] Effect of -pc_gamg_threshold vs PETSc version

[Jeremy Theler]

Hi Mark. So glad you answered.

> 0) what is your test problem? eg, 3D Lapacian with Q1 finite
> elements.

I said in my first email it was linear elasticty (and I gave a link
where you can see the geometry, BCs, etc.) but I did not specifty
further details.
It is linear elasticity with displacement-based FEM formulation using
unstructured curved 10-noded tetrahedra.

The matrix is marked as SPD with MatSetOption() and the solver is
indeed CG and not the default GMRES.

> First, you can get GAMG diagnostics by running with '-info :pc' and
> grep on GAMG.

Great advice. Now I have a lot more of information but I'm not sure how
to analyze it. Find attached for each combination of threshold and
PETSc version the output of -info :pc -ksp_monitor -ksp_view

In general it looks like 3.18 and 3.19 have less KSP iterations than
3.17 but the overall time is larger.

> Anyway, sorry for the changes. 
> I hate changing GAMG for this reason and I hate AMG for this reason!

No need to apologize, I just want to better understand how to better
exploit your code!

Thanks
--
jeremy

> 
> Thanks,
> Mark
> 
> 
> 
> On Thu, Apr 13, 2023 at 8:17 AM Jeremy Theler 
> wrote:
> > When using GAMG+cg for linear elasticity and providing the near
> > nullspace computed by MatNullSpaceCreateRigidBody(), I used to find
> > "experimentally" that a small value of -pc_gamg_threshold in the
> > order
> > of 0.0001 would slightly decrease the solve time.
> > 
> > Starting with 3.18, I started seeing that any positive value for
> > the
> > treshold would increase the solve time. I did a quick parametric
> > (serial) run solving an elastic problem with a matrix size of
> > approx
> > 570k x 570k for different values of GAMG threshold and different
> > PETSc
> > versions (compiled with the same compiler, options and flags).
> > 
> > I noted that
> > 
> >  1. starting from 3.18, a threshold of 0.0001 that used to improve
> > the
> > speed now worsens it. 
> >  2. PETSc 3.17 looks like a "sweet spot" of speed
> > 
> > I would like to hear any comments you might have.
> > 
> > The wall time shown includes the time needed to read the mesh and
> > assemble the stiffness matrix. It is a refined version of the
> > NAFEMS
> > LE10 benchmark described here:
> > https://seamplex.com/feenox/examples/mechanical.html#nafems-le10-
> > thick-plate-pressure-benchmark
> > 
> > If you want, I could dump the matrix, rhs and near nullspace
> > vectors
> > and share them.
> > 
> > --
> > jeremy theler
> > 



log.tar.gz
Description: application/compressed-tar

[petsc-users] Effect of -pc_gamg_threshold vs PETSc version

[Jeremy Theler]

When using GAMG+cg for linear elasticity and providing the near
nullspace computed by MatNullSpaceCreateRigidBody(), I used to find
"experimentally" that a small value of -pc_gamg_threshold in the order
of 0.0001 would slightly decrease the solve time.

Starting with 3.18, I started seeing that any positive value for the
treshold would increase the solve time. I did a quick parametric
(serial) run solving an elastic problem with a matrix size of approx
570k x 570k for different values of GAMG threshold and different PETSc
versions (compiled with the same compiler, options and flags).

I noted that

 1. starting from 3.18, a threshold of 0.0001 that used to improve the
speed now worsens it. 
 2. PETSc 3.17 looks like a "sweet spot" of speed

I would like to hear any comments you might have.

The wall time shown includes the time needed to read the mesh and
assemble the stiffness matrix. It is a refined version of the NAFEMS
LE10 benchmark described here:
https://seamplex.com/feenox/examples/mechanical.html#nafems-le10-thick-plate-pressure-benchmark

If you want, I could dump the matrix, rhs and near nullspace vectors
and share them.

--
jeremy theler



threshold.pdf
Description: Adobe PDF document
15 30.06
16 30.44
17 26.29
18 28.61
19 29.65
15 29.23
16 29.54
17 24.70
18 29.44
19 30.58
15 30.11
16 30.51
17 25.80
18 32.68
19 33.78
15 31.68
16 31.96
17 26.98
18 43.36
19 44.24
15 36.74
16 37.06
17 31.96
18 69.54
19 70.14

[petsc-users] API cal to set mg_levels_pc_type

[Jeremy Theler]

The default smoothing PC changed from sor to jacobi in 3.17. The
Changelog says the old behavior can be recovered by using
-mg_levels_pc_type sor.

1. Is there a way to set the mg_levels_pc_type via an API call?
2. Are there any changes in efficiency expected with this new PC?


Regards
--
jeremy theler

Re: [petsc-users] set MUMPS OOC_TMPDIR

[Jeremy Theler]

I've known at least one other person that had this very same confusion.
Maybe it's worth adding Barry's explanation to the manual page:

This function can be called BEFORE PetscInitialize(), but it does not
NEED to be called before PetscInitialize(). You should be able to call
PetscOptionsSetValue() anytime you want.




On Mon, 2022-02-28 at 16:39 -0800, Sam Guo wrote:
> Yes, " This function can be called BEFORE PetscInitialize()" confused
> me. Thanks for the clarification. 
> 
> On Mon, Feb 28, 2022 at 4:38 PM Barry Smith  wrote:
> > 
> >   You should be able to call PetscOptionsSetValue() anytime you
> > want, as I said between different uses of MUMPS you can call it to
> > use different directories.
> > 
> >   Perhaps this confused you?
> > 
> >      Note:
> >    This function can be called BEFORE PetscInitialize()
> > 
> >   It is one of the very few functions that can be called before
> > PetscInitialize() but it does not NEED to be called before
> > PetscInitialize().
> > 
> >   Barry

Re: [petsc-users] HDF5 timestepping in PETSc 3.16

[Jeremy Theler]

On Thu, 2021-10-21 at 13:04 -0400, Matthew Knepley wrote:
> On Tue, Oct 19, 2021 at 6:12 AM Matthew Knepley 
> wrote:
> > On Mon, Oct 18, 2021 at 10:35 PM Adrian Croucher
> >  wrote:
> > > Any response on this?
> > > 
> > > This is a bit of a showstopper for me - I can't upgrade to PETSc
> > > 3.16 if 
> > > it does not allow my users to read their HDF5 files created using
> > > earlier versions of PETSc.
> > > 
> > > So far I can't see a workaround. Possibly the timestepping
> > > functions 
> > > need some kind of optional parameter to specify what the default 
> > > timestepping attribute should be, if it's not present in the file
> > 
> > I think you are right. We should always write the attribute, but
> > have it be false. We should
> > interpret a missing attribute as an old file.
> > 
> Okay, I think I have it. Can you look at this branch?
> 
>   https://gitlab.com/petsc/petsc/-/merge_requests/4483
> 
> There is now an option that lets you set the default timestepping
> behavior
> 
>   -viewer_hdf5_default_timestepping
> 
> I think that is what you want.

I'd like to rely on PetscViewerHDF5SetDefaultTimestepping() to provide
backwards compatibility as well. This branch has been merged into
master back in November but never made it to the stable v3.16 releases.

Can you guys check please?


Thanks
--
jeremy theler
www.seamplex.com

Re: [petsc-users] DM misuse causes massive memory leak?

[Jeremy Theler]

Hola Jesús

You were lucky you got 20x. See

https://petsc.org/release/docs/manualpages/Mat/MatSeqAIJSetPreallocation.html
https://petsc.org/release/docs/manualpages/Mat/MatMPIAIJSetPreallocation.html

Quote:

"For large problems you MUST preallocate memory or you will get
TERRIBLE performance, see the users' manual chapter on matrices. "

--
jeremy


On Wed, 2021-12-29 at 22:12 +, Ferrand, Jesus A. wrote:
> Dear PETSc Team:
> 
> I have a question about DM andPetscSection. Say I import a mesh (for
> FEM purposes) and create aDMPlex for it. I then usePetscSections to
> set degrees of freedom per "point" (by point I mean vertices, lines,
> faces, and cells). I then usePetscSectionGetStorageSize() to get the
> size of the global stiffness matrix (K) needed for my FEM problem.
> One last detail, this K I populate inside a rather large loop using
> an element stiffness matrix function of my own. Instead of
> usingDMCreateMatrix(), I manually created aMat using MatCreate(),
> MatSetType(),MatSetSizes(), and MatSetUp(). I come to find that said
> loop is painfully slow when I use the manually created matrix, but
> 20x faster when I use the Mat coming out ofDMCreateMatrix().
> 
> My question is then: Is the manual Mat a noob mistake and is it
> somehow creating a memory leak with K? Just in case it's something
> else I'm attaching the code. The loop that populates K is between
> lines 221 and 278. Anything related to DM, DMPlex, and PetscSection
> is between lines 117 and 180.
> 
> Machine Type: HP Laptop
> C-compiler: Gnu C
> OS: Ubuntu 20.04
> PETSc version: 3.16.0
> MPI Implementation: MPICH
> 
> Hope you all had a Merry Christmas and that you will have a happy and
> productive New Year. :D
> 
> Sincerely:
> J.A. Ferrand
> Embry-Riddle Aeronautical University - Daytona Beach FL
> M.Sc. Aerospace Engineering | May 2022
> B.Sc. Aerospace Engineering
> B.Sc. Computational Mathematics
>  
> Sigma Gamma Tau 
> Tau Beta Pi
> Honors Program
>  
> Phone: (386)-843-1829
> Email(s): ferra...@my.erau.edu
>    jesus.ferr...@gmail.com

Re: [petsc-users] How to read/write a HDF5 file using petsc4py ?

[Jeremy Theler]

On Fri, 2021-12-10 at 10:02 +0100, Quentin Chevalier wrote:
> root@container:/home/shared# echo $PETSCH_ARCH $PETSCH_DIR

PETSCH_ARCH -> PETSC_ARCH
PETSCH_DIR -> PETSC_DIR

Re: [petsc-users] Finding which cell an arbitrary point belongs to in DMPlex

[Jeremy Theler]

On Thu, 2020-09-24 at 18:30 -0400, Matthew Knepley wrote:
> > There is also DMPlexFindVertices() which finds the nearest vertex
> > to
> > > the given coords in the given radius.
> > 
> > At first I had understood that this function performed a nearest-
> > neighbor search but after a closer look it sweeps the local DM and
> > marks whether the sought points coincide with a mesh node within
> > eps or
> > not. Neat.

This DMPlexFindVertices() sweeps over DMGetCoordinatesLocal() which
returns both the local and ghost coordinates, so at the end of the day
I might get more than one process claiming to have found the same node.

How can I ignore ghost points so each vertex actually belongs to the
process that found it?


> > > I wrote it some time ago mainly for debug purposes. It uses just
> > > brute force. I'm not sure it deserves to exist :-) Maybe we
> > should
> > > somehow merge these functionalities.
> > 
> > It works, although a kd-tree-based search would be far more
> > efficient
> > than a full sweep over the DM.
> 
> We should not need to do that. LocatePoints() does not sweep the
> mesh.
> It just does grid hashing. kd is a little better with really
> irregular distributions,
> but hashing should be fine.

Yes, it seems to be pretty efficent (although there is no support for
3D so far).

Two more things about plexgeometry.c:

 1. shouldn't line 224 in DMPlexLocatePoint_Simplex_3D_Internal()
compare against -eps instead of against zero as donde in line 145 in
DMPlexLocatePoint_Simplex_2D_Internal()?

 2. wouldn't it be better to replace DMGetDimension() by
DMGetCoordinateDim() in line 45 inside DMPlexFindVertices? I have a 2D
mesh with 3D coordinates and PetscUnlikely() is triggered.

Thanks
--
jeremy

Re: [petsc-users] Finding which cell an arbitrary point belongs to in DMPlex

[Jeremy Theler]

On Wed, 2020-09-16 at 14:29 +, Hapla  Vaclav wrote:
> There is also DMPlexFindVertices() which finds the nearest vertex to
> the given coords in the given radius.

At first I had understood that this function performed a nearest-
neighbor search but after a closer look it sweeps the local DM and
marks whether the sought points coincide with a mesh node within eps or
not. Neat.


> You can then get support or its transitive closure for that vertex.

Not directly because in general the sought points will not coincide
with a mesh node, but a combination of this function and
DMLocatePoints() seems to do the trick.  

> I wrote it some time ago mainly for debug purposes. It uses just
> brute force. I'm not sure it deserves to exist :-) Maybe we should
> somehow merge these functionalities.

It works, although a kd-tree-based search would be far more efficient
than a full sweep over the DM.

Thanks
--
jeremy

> 
> Thanks,
> 
> Vaclav
> 
> > On 16 Sep 2020, at 01:44, Matthew Knepley 
> > wrote:
> > 
> > On Tue, Sep 15, 2020 at 6:18 PM Jeremy Theler 
> > wrote:
> > > On Mon, 2020-09-14 at 20:28 -0400, Matthew Knepley wrote:
> > > > On Mon, Sep 14, 2020 at 6:15 PM Jeremy Theler <
> > > jer...@seamplex.com>
> > > > wrote:
> > > > > Hello all
> > > > > 
> > > > > Say I have a fully-interpolated 3D DMPlex and a point with
> > > > > arbitrary
> > > > > coordinates x,y,z. What's the most efficient way to know
> > > which cell
> > > > > this point belongs to in parallel? Cells can be either tets
> > > or
> > > > > hexes.
> > > > 
> > > > I should make a tutorial on this, but have not had time so far.
> > > 
> > > Thank you very much for this mini-tutorial.
> > > 
> > > > 
> > > > The intention is that you use
> > > > 
> > > >   
> > > > 
> > > https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DM/DMLocatePoints.html
> > > > 
> > > > This will just brute force search unless you also give
> > > > 
> > > >   -dm_plex_hash_location
> > > 
> > > Well, for a 3D DMplex PETSc (and git blame) tells me that you
> > > "have
> > > only coded this for 2D." :-)
> > > 
> > 
> > Crap. I need to do 3D. It's not hard, just work.
> >  
> > > > which builds a grid hash to accelerate it. I should probably
> > > expose
> > > > 
> > > >   DMPlexLocatePoint_Internal()
> > > > 
> > > > which handles the single cell queries. If you just had one
> > > point,
> > > > that might make it simpler,
> > > > although you would still write your own loop.
> > > 
> > > I see that DMLocatePoints() loops over all the cells until it
> > > finds the
> > > right one. I was thinking about finding first the nearest vertex
> > > to the
> > > point and then sweeping over all the cells that share this vertex
> > > testing for DMPlexLocatePoint_Internal(). The nearest node ought
> > > to be
> > > found using an octree or similar. Any direction regarding this
> > > idea?
> > > 
> > 
> > So you can imagine both a topological search and a geometric
> > search. Generally, people want geometric.
> > The geometric hash we use is just to bin elements on a regular
> > grid.
> >  
> > > >  If your intention is to interpolate a field at these
> > > > locations, I created
> > > > 
> > > >   
> > > > 
> > > https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/SNES/DMInterpolationCreate.html
> > > > 
> > > > which no one but me uses so far, but I think it is convenient.
> > > 
> > > Any other example apart from src/snes/tutorials/ex63.c?
> > > 
> > 
> > That is the only one in PETSc. The PyLith code uses this to
> > interpolate to seismic stations.
> > 
> >   Thanks,
> > 
> >  Matt
> >  
> > > Thank you.
> > > 
> > > > 
> > > >   Thanks,
> > > > 
> > > > Matt
> > > >  
> > > > > Regards
> > > > > --
> > > > > jeremy theler
> > > > > www.seamplex.com
> > > > > 
> > > > > 
> > > > 
> > > > 
> > > > 
> > > 
> > > 
> > 
> > 
> > -- 
> > 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] Finding which cell an arbitrary point belongs to in DMPlex

[Jeremy Theler]

On Mon, 2020-09-14 at 20:28 -0400, Matthew Knepley wrote:
> On Mon, Sep 14, 2020 at 6:15 PM Jeremy Theler 
> wrote:
> > Hello all
> > 
> > Say I have a fully-interpolated 3D DMPlex and a point with
> > arbitrary
> > coordinates x,y,z. What's the most efficient way to know which cell
> > this point belongs to in parallel? Cells can be either tets or
> > hexes.
> 
> I should make a tutorial on this, but have not had time so far.

Thank you very much for this mini-tutorial.

> 
> The intention is that you use
> 
>   
> https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DM/DMLocatePoints.html
> 
> This will just brute force search unless you also give
> 
>   -dm_plex_hash_location

Well, for a 3D DMplex PETSc (and git blame) tells me that you "have
only coded this for 2D." :-)


> which builds a grid hash to accelerate it. I should probably expose
> 
>   DMPlexLocatePoint_Internal()
> 
> which handles the single cell queries. If you just had one point,
> that might make it simpler,
> although you would still write your own loop.

I see that DMLocatePoints() loops over all the cells until it finds the
right one. I was thinking about finding first the nearest vertex to the
point and then sweeping over all the cells that share this vertex
testing for DMPlexLocatePoint_Internal(). The nearest node ought to be
found using an octree or similar. Any direction regarding this idea?


>  If your intention is to interpolate a field at these
> locations, I created
> 
>   
> https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/SNES/DMInterpolationCreate.html
> 
> which no one but me uses so far, but I think it is convenient.

Any other example apart from src/snes/tutorials/ex63.c?


Thank you.

> 
>   Thanks,
> 
> Matt
>  
> > Regards
> > --
> > jeremy theler
> > www.seamplex.com
> > 
> > 
> 
> 
> 

[petsc-users] Finding which cell an arbitrary point belongs to in DMPlex

[Jeremy Theler]

Hello all

Say I have a fully-interpolated 3D DMPlex and a point with arbitrary
coordinates x,y,z. What's the most efficient way to know which cell
this point belongs to in parallel? Cells can be either tets or hexes.

Regards
--
jeremy theler
www.seamplex.com

Re: [petsc-users] BCs for a EPS solver

[Jeremy Theler]

The usual trick is to set ones in one matrix and zeros in the other
one.


On Mon, 2020-02-17 at 12:35 -0600, Emmanuel Ayala wrote:
> Hi everyone,
> 
> I have an eigenvalue problem where I need to apply BCs to the
> stiffness and mass matrix. 
> 
> Usually, for KSP solver, it is enough to set to zero the rows and
> columns related to the boundary conditions. I used to apply it with
> MatZeroRowsColumns, with a 1s on the diagonal. Then the solver works
> well.
> 
> There is something similar to KSP for EPS solver ? 
> 
> I already used MatZeroRowsColumns (for EPS solver), with a 1s on the
> diagonal, and I got wrong result.
> 
> Kind regards.
> 
> 
> 
> 

Re: [petsc-users] Product of matrix row times a vector

[Jeremy Theler]

On Thu, 2020-01-30 at 18:05 -0500, Matthew Knepley wrote:
> On Thu, Jan 30, 2020 at 6:04 PM Jeremy Theler 
> wrote:
> > On Thu, 2020-01-30 at 21:10 +, Smith, Barry F. wrote:
> > 
> > >   MatGetSubMatrix() and then do the product on the sub matrix
> > then
> > > VecSum
> > > 
> > 
> > Ok, I have it working in a single-processor and throws the expected
> > value. Yet I have a segfault in parallel when I ask for the IS
> > corresponding to the rows. When I call this instruction in parallel
> > I
> > get a segfault (I can post the full debug output if needed):
> > 
> > ISCreateStride(PETSC_COMM_WORLD, size_local, first_row, 1,
> > _cols);
> > 
> > I tried also
> > 
> > ISCreateStride(PETSC_COMM_WORLD, size_global, 0, 1, _cols));
> > 
> > but it also fails with the same segfault.
> > 
> > What am I getting wrong?
> 
> Show the entire error, including stack.
> 

I run in two processes with start_on_debugger. Main terminal says

malloc(): corrupted top size


and the gdb output with the stack trace is


Attaching to program: /home/gtheler/codigos/wasora-suite/fino/examples/fino, 
process 31192
[New LWP 31196]
[New LWP 31198]
[Thread debugging using libthread_db enabled]
Using host libthread_db library "/lib/x86_64-linux-gnu/libthread_db.so.1".
0x7fd236019720 in __GI___nanosleep (
requested_time=requested_time@entry=0x7ffdf9838190, 
remaining=remaining@entry=0x7ffdf9838190)
at ../sysdeps/unix/sysv/linux/nanosleep.c:28
28  ../sysdeps/unix/sysv/linux/nanosleep.c: No such file or directory.
(gdb) c
Continuing.
[New Thread 0x7fd233034700 (LWP 31212)]
Thread 1 "fino" received signal SIGABRT, Aborted.
__GI_raise (sig=sig@entry=6) at ../sysdeps/unix/sysv/linux/raise.c:50
50  ../sysdeps/unix/sysv/linux/raise.c: No such file or directory.
(gdb) where
#0  __GI_raise (sig=sig@entry=6) at ../sysdeps/unix/sysv/linux/raise.c:50
#1  0x7fd235f75535 in __GI_abort () at abort.c:79
#2  0x7fd235fcc508 in __libc_message (action=action@entry=do_abort, 
fmt=fmt@entry=0x7fd2360d728d "%s\n") at ../sysdeps/posix/libc_fatal.c:181
#3  0x7fd235fd2c1a in malloc_printerr (
str=str@entry=0x7fd2360d5518 "malloc(): corrupted top size")
at malloc.c:5341
#4  0x7fd235fd620d in _int_malloc (
av=av@entry=0x7fd23610ec40 , bytes=bytes@entry=2372)
at malloc.c:4099
#5  0x7fd235fd756a in __GI___libc_malloc (bytes=2372) at malloc.c:3057
#6  0x7fd2377dd42d in PetscMallocAlign (mem=2372, clear=PETSC_TRUE, 
line=37, func=0x7fd239450408 <__func__.14861> "ISCreate", 
file=0x7fd239450298 
"/home/gtheler/libs/petsc-3.12.3/src/vec/is/is/interface/isreg.c", 
result=0x7ffdf983c618)
at /home/gtheler/libs/petsc-3.12.3/src/sys/memory/mal.c:49
#7  0x7fd2377e0875 in PetscTrMallocDefault (a=768, clear=PETSC_TRUE, 
lineno=37, function=0x7fd239450408 <__func__.14861> "ISCreate", 
filename=0x7fd239450298 
"/home/gtheler/libs/petsc-3.12.3/src/vec/is/is/interface/isreg.c", 
result=0x7ffdf983c8b0)
at /home/gtheler/libs/petsc-3.12.3/src/sys/memory/mtr.c:164
#8  0x7fd2377dedf0 in PetscMallocA (n=1, clear=PETSC_TRUE, lineno=37, 
function=0x7fd239450408 <__func__.14861> "ISCreate", 
--Type  for more, q to quit, c to continue without paging--
filename=0x7fd239450298 
"/home/gtheler/libs/petsc-3.12.3/src/vec/is/is/interface/isreg.c", bytes0=768, 
ptr0=0x7ffdf983c8b0)
at /home/gtheler/libs/petsc-3.12.3/src/sys/memory/mal.c:422
#9  0x7fd237c67760 in ISCreate (comm=0x556dbe8c9ee0 , 
is=0x7ffdf983c8b0)
at /home/gtheler/libs/petsc-3.12.3/src/vec/is/is/interface/isreg.c:37
#10 0x7fd237c4589a in ISCreateStride (
comm=0x556dbe8c9ee0 , n=70188, first=70191, step=1, 
is=0x7ffdf983c8b0)
#11 0x556dbe82aab4 in fino_instruction_reaction (arg=0x556dbf8b4460)
at ./reactions.c:72
#12 0x556dbe8a8aa2 in wasora_step (whence=0) at ../wasora/src/wasora.c:406
#13 0x556dbe8a81b3 in wasora_standard_run () at ../wasora/src/wasora.c:206
#14 0x556dbe8a80a5 in main (argc=3, argv=0x7ffdf983cbd8)
at ../wasora/src/wasora.c:166
(gdb) 

Re: [petsc-users] Product of matrix row times a vector

[Jeremy Theler]

On Thu, 2020-01-30 at 21:10 +, Smith, Barry F. wrote:

>   MatGetSubMatrix() and then do the product on the sub matrix then
> VecSum
> 

Ok, I have it working in a single-processor and throws the expected
value. Yet I have a segfault in parallel when I ask for the IS
corresponding to the rows. When I call this instruction in parallel I
get a segfault (I can post the full debug output if needed):

ISCreateStride(PETSC_COMM_WORLD, size_local, first_row, 1, _cols);

I tried also

ISCreateStride(PETSC_COMM_WORLD, size_global, 0, 1, _cols));

but it also fails with the same segfault.

What am I getting wrong?

--
jeremy

[petsc-users] Product of matrix row times a vector

[Jeremy Theler]

Sorry if this is basic, but I cannot figure out how to do it in
parallel and I'd rather not say how I do it in single-processor mode
because I would be ashamed.

Say I have a matrix and I want to multiply a row times a vector to
obtain a scalar. Actually I would like to choose some rows, multiply
each of them by the vector and then add the scalars up. Or, conversely,
sum all the rows colmunwise and then multiply the sum by a vector.

How can I do this?


Regards
--
jeremy theler
www.seamplex.com

[petsc-users] Internal product through a matrix norm

[Jeremy Theler]

Sorry for the basic question, but here it goes.
Say I have a vector u and a matrix K and I want to compute the scalar

e = u^T K u

(for example the strain energy if u are displacements a K is the
stiffness matrix).

Is there anything better (both in elegance and efficiency) than doing
this?

PetscScalar e;
Vec Kx;
  
VecDuplicate(x, );
MatMult(K, x, Kx);
VecDot(x, Kx, );


--
jeremy theler
www.seamplex.com

Re: [petsc-users] error handling

[Jeremy Theler]

Dear Sam

Probably you are already aware of the following paragraph, but just in
case. Quote from 
https://www.gnu.org/prep/standards/standards.html#Memory-Usage


Memory analysis tools such as valgrind can be useful, but don’t
complicate a program merely to avoid their false alarms. For example,
if memory is used until just before a process exits, don’t free it
simply to silence such a tool. 


Regards
--
jeremy theler
www.seamplex.com



On Tue, 2020-01-21 at 08:49 -0800, Sam Guo wrote:
> I use PETSc from my application. Sounds you are saying I just treat
> ierr!=0 as an system error and no need to call Destroy functions.
> 
> On Tuesday, January 21, 2020, Smith, Barry F. 
> wrote:
> > 
> > > On Jan 20, 2020, at 6:32 PM, Sam Guo 
> > wrote:
> > > 
> > > Hi Barry,
> > >   I understand ierr != 0 means  something catastrophic. I just
> > want to release all memory before I exit PETSc.
> > 
> >In general not possible. If you run with the debug version and
> > -malloc_debug it is possible but because of the unknown error it
> > could be that the releasing of the memory causes a real crash.
> > 
> >Is your main concern when you use PETSc for a large problem and
> > it errors because it is "out of memory"?
> > 
> >Barry
> > 
> > 
> > > 
> > > Thanks,
> > > Sam
> > > 
> > > On Mon, Jan 20, 2020 at 4:06 PM Smith, Barry F. <
> > bsm...@mcs.anl.gov> wrote:
> > > 
> > >   Sam,
> > > 
> > > I am not sure what your goal is but PETSc error return codes
> > are error return codes not exceptions. They mean that something
> > catastrophic happened and there is no recovery. 
> > > 
> > > Note that PETSc solvers do not return nonzero error codes on
> > failure to converge etc. You call, for example,
> > KPSGetConvergedReason() after a KSP solve to see if it has failed,
> > this is not a catastrophic failure. If a MatCreate() or any other
> > call returns a nonzero ierr the game is up, you cannot continue
> > running PETSc.
> > > 
> > >Barry
> > > 
> > > 
> > > > On Jan 20, 2020, at 5:41 PM, Matthew Knepley  > > wrote:
> > > > 
> > > > Not if you initialize the pointers to zero: Mat A = NULL.
> > > > 
> > > >Matt
> > > > 
> > > > On Mon, Jan 20, 2020 at 6:31 PM Sam Guo 
> > wrote:
> > > > I mean MatDestroy.
> > > > 
> > > > On Mon, Jan 20, 2020 at 3:28 PM Sam Guo 
> > wrote:
> > > > Does it hurt to call Destroy function without calling
> > CreateFunction? For example 
> > > > Mat A, B;
> > > > PetscErrorCode  ierr1, ierr2;
> > > > ierr1 = MatCreate(PETSC_COMM_WORLD,);
> > > > if(ierr1 == 0)
> > > > {
> > > >   ierr2 = MatCreate(PETSC_COMM_WORLD
> > > > ,);
> > > > 
> > > > }
> > > > if(ierr1 !=0 || ierr2 != 0)
> > > > {
> > > >   Destroy();
> > > >   Destroy(); // if ierr1 !=0, MatCreat is not called on B.
> > Does it hurt to call Destroy B here?
> > > > }
> > > > 
> > > > 
> > > > 
> > > > On Mon, Jan 20, 2020 at 11:11 AM Dave May <
> > dave.mayhe...@gmail.com> wrote:
> > > > 
> > > > 
> > > > On Mon 20. Jan 2020 at 19:47, Sam Guo 
> > wrote:
> > > > Can I assume if there is MatCreat or VecCreate, I should clean
> > up the memory myself?
> > > > 
> > > > Yes. You will need to call the matching Destroy function.
> > > > 
> > > > 
> > > > 
> > > > On Mon, Jan 20, 2020 at 10:45 AM Sam Guo  > > wrote:
> > > > I only include the first few lines of SLEPc example. What about
> > following
> > > >   ierr = MatCreate(PETSC_COMM_WORLD,);CHKERRQ(ierr); 
> > > >   ierr =
> > MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n);CHKERRQ(ierr);  
> > > > Is there any memory  lost?
> > > > 
> > > > On Mon, Jan 20, 2020 at 10:41 AM Dave May <
> > dave.mayhe...@gmail.com> wrote:
> > > > 
> > > > 
> > > > On Mon 20. Jan 2020 at 19:39, Sam Guo 
> > wrote:
> > > > I don't have a specific case yet. Currently every call of PETSc
> > is checked. If ierr is not zero, print the error and return. For
> > example,
> > > >Mat A; /* problem matrix */ 
> > > >EPS eps; /* eigenproblem solver c

Re: [petsc-users] pc_gamg_threshol

[Jeremy Theler]

Yes, I read that page and it was that paragraph that made me want to
learn more.

For example, that pages says:

“-pc_gamg_threshold 0.0 is the most robust option (the reason for this
is not obvious) ...”


Where can I find more math-based background on this subject? I mean,
some text that describes the methods and not just the implementation as
the source code at gamg/util.c so I can better understand what is going
on.


Thanks



-- 
Jeremy Theler
www.seamplex.com



On Thu, 2017-01-05 at 09:18 -0500, Mark Adams wrote:
> You want the bottom of page 84 in the manual.
> 
> On Wed, Jan 4, 2017 at 4:33 PM, Barry Smith <bsm...@mcs.anl.gov>
> wrote:
> 
>The manual page gives a high-level description
> 
> http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/PC/PCGAMGSetThreshold.html
>  the exact details can be found in the code here 
> http://www.mcs.anl.gov/petsc/petsc-dev/src/ksp/pc/impls/gamg/util.c.html#PCGAMGFilterGraph
>   I'm adding a link from the former to the later in the documentation.
> 
>Barry
> 
> 
>     
> > On Jan 4, 2017, at 3:16 PM, Jeremy Theler
> <jer...@seamplex.com> wrote:
> >
> > * Any reference to what pc_gamg_treshold means and/or does?
> >
> >
> >
> > On Wed, 2017-01-04 at 18:13 -0300, Jeremy Theler wrote:
> >> Hi! Any reference to what does -pc_gamg_threshold mean
> and/or?
> >>
> >
> 
> 
> 
> 

Re: [petsc-users] pc_gamg_threshol

[Jeremy Theler]

* Any reference to what pc_gamg_treshold means and/or does?



On Wed, 2017-01-04 at 18:13 -0300, Jeremy Theler wrote:
> Hi! Any reference to what does -pc_gamg_threshold mean and/or?
> 

[petsc-users] pc_gamg_threshol

[Jeremy Theler]

Hi! Any reference to what does -pc_gamg_threshold mean and/or?

-- 
Jeremy Theler
www.seamplex.com

[petsc-users] getting the near nullspace from PCSetCoordinates

[Jeremy Theler]

Hi all

I want to check that the near nullspace I provide to GAMG gives "almost
null vectors" when multiplying each vector in the near nullspace against
the matrix problem.

This way I can check that the unknown ordering I am using is consistent,
for example using by MatNullSpaceCreateRigidBody() or by computing the
nullspace by myself.


The thing is I do not know how I can get the nullspace object after
calling PCSetCoordinates(). It gets a pointer to the PC object, but
MatGetNearNullSpace() needs the matrix object. I assume at some point
the matrix and the PC get linked, but when I ask
MatGetNearNullSpace(matrix) passing the problem matrix after setting
PCSetCoordinates(pc) I get:

error: PETSc error 85-0 'Null Object: Parameter # 1'
in /home/gtheler/libs/petsc-3.7.4/src/mat/interface/matnull.c
MatNullSpaceGetVecs:64


thanks

-- 
Jeremy Theler
www.seamplex.com

Re: [petsc-users] GAMG

[Jeremy Theler]

On Mon, 2016-10-31 at 08:44 -0600, Jed Brown wrote:

> > After understanding Matt's point about the near nullspace (and reading
> > some interesting comments from Jed on scicomp stackexchange) I did built
> > my own vectors (I had to take a look at MatNullSpaceCreateRigidBody()
> > because I found out by running the code the nullspace should be an
> > orthonormal basis, it should say so in the docs).
> 
> Where?
> "vecs - the vectors that span the null space (excluding the constant 
> vector); these vectors must be orthonormal."
> https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatNullSpaceCreate.html

ok, I might have passed on that but I started with 

http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatSetNearNullSpace.html

that says “attaches a null space to a matrix, which is often the null
space (rigid body modes) of the operator without boundary conditions
This null space will be used to provide near null space vectors to a
multigrid preconditioner built from this matrix.”

It wouldn't hurt to remind dumb users like me that “...it is often the
set of _orthonormalized_ rigid body modes...”


> And if you run in debug mode (default), as you always should until you
> are confident that your code is correct, MatNullSpaceCreate tests that
> your vectors are orthonormal.

That's how I realized I needed to normalize. Then I found
MatNullSpaceCreateRigidBody() and copied the code to orthogonalize.

Wouldn't it be better to orthonormalize inside MatSetNullSpace()? I bet
an orthogonalization from PETSc's code would beat any user-side code.

> > Now, there are some results I do not understand. I tried these six
> > combinations:
> >
> > order  near-nullspace   iterationsnorm
> > -  --   --
> > unknownexplicit 101.6e-6
> > unknownPCSetCoordinates 151.7e-7
> > unknownnone 152.4e-7
> > node   explicit fails with error -11
> > node   PCSetCoordinates fails with error -11
> > node   none 133.8e-7
> 
> Did you set a block size for the "node-based" orderings?  Are you sure
> the above is labeled correctly?  Anyway, PCSetCoordinates uses
> "node-based" ordering.  Implementation performance will generally be
> better with node-based ordering -- it has better memory streaming and
> cache behavior.

Yes. Indeed, when I save the stiffnes matrix as a binary file I get
a .info file that contains

-matload_block_size 3

The labeling is right, I re-checked. That's the funny part, I can't get
GAMG to work with PCSetCoordinates (which BTW, I think its documentation
does not address the issue of DOF ordering).

Any idea of what can be happening to me?


> The AIJ matrix format will also automatically do an "inode" optimization
> to reduce memory bandwidth and enable block smoothing (default
> configuration uses SOR smoothing).  You can use -mat_no_inode to try
> turning that off.


That option does not make any difference.

> 
> > Error -11 is 
> > PETSc's linear solver did not converge with reason
> > 'DIVERGED_PCSETUP_FAILED' (-11)
> Isn't there an actual error message?

Sorry, KSPGetConvergedReason() returns -11 and then my code prints that
error string. Find attached the output with -info.

Thanks
--
jeremy



[0] PetscInitialize(): PETSc successfully started: number of processors = 1
[0] PetscGetHostName(): Rejecting domainname, likely is NIS tom.(none)
[0] PetscInitialize(): Running on machine: tom
[0] SlepcInitialize(): SLEPc successfully started
697 3611
[0] PetscCommDuplicate(): Duplicating a communicator 2 2 max tags = 1
[0] MatAssemblyEnd_SeqAIJ(): Matrix size: 2091 X 2091; storage space: 835506 
unneeded,74079 used
[0] MatAssemblyEnd_SeqAIJ(): Number of mallocs during MatSetValues() is 0
[0] MatAssemblyEnd_SeqAIJ(): Maximum nonzeros in any row is 81
[0] MatCheckCompressedRow(): Found the ratio (num_zerorows 0)/(num_localrows 
2091) < 0.6. Do not use CompressedRow routines.
[0] MatSeqAIJCheckInode(): Found 697 nodes of 2091. Limit used: 5. Using Inode 
routines
[0] MatAssemblyEnd_SeqAIJ(): Matrix size: 2091 X 2091; storage space: 0 
unneeded,74079 used
[0] MatAssemblyEnd_SeqAIJ(): Number of mallocs during MatSetValues() is 0
[0] MatAssemblyEnd_SeqAIJ(): Maximum nonzeros in any row is 81
[0] MatCheckCompressedRow(): Found the ratio (num_zerorows 0)/(num_localrows 
2091) < 0.6. Do not use CompressedRow routines.
[0] MatAssemblyEnd_SeqAIJ(): Matrix size: 2091 X 2091; storage space: 0 
unneeded,74079 used
[0] MatAssemblyEnd_SeqAIJ(): Number of mallocs during MatSetValues() is 0
[0] MatAssemblyEnd_SeqAIJ(): Maximum nonzeros in any row is 81
[0] MatCheckCompressedRow(): Found the ratio (num_zerorows 0)/(num_localrows 
2091) < 0.6. Do not use CompressedRow routines.
[0] PCSetUp(): Setting up PC for first time
[0] PCSetUp_GAMG(): level 0) N=2091, n data rows=3, n data cols=6, nnz/row 

Re: [petsc-users] GAMG

[jeremy theler]

I will try these options in a couple of hours (I have to go out now). I
forgot to mention that the geometry has revolution symmetry around the z
axis (just the geometry, not the problem because it has a non-symmetric
temperature distribution).
I am solving with only one proc, there are approx 50k nodes so 150k dofs.
Thanks again.

On Fri, Oct 28, 2016, 11:07 Mark Adams <mfad...@lbl.gov> wrote:

> Also, try solving the problem with a one level iterative method and
> Chebyshev, like:
>
> -ksp_type chebyshev
> -pc_type jacobi
>
> It will take a long time to solve but I just want to see if it has the
> same error.
>
>
> On Fri, Oct 28, 2016 at 10:04 AM, Mark Adams <mfad...@lbl.gov> wrote:
>
> GAMG's eigen estimator worked but the values are very high.  You have very
> low number of equations per processor, is this a thin body? Are the
> elements badly stretched?
>
> Do this again with these parameters:
>
> -mg_levels_ksp_type chebyshev
> -mg_levels_esteig_ksp_type cg
> -mg_levels_esteig_ksp_max_it 10
> ​​
> -mg_levels_ksp_chebyshev_esteig 0,.1,0,1.05
> -gamg_est_ksp_type cg
>
>
> On Fri, Oct 28, 2016 at 9:48 AM, Jeremy Theler <jer...@seamplex.com>
> wrote:
>
> On Fri, 2016-10-28 at 09:46 -0400, Mark Adams wrote:
> > Please run with -info and grep on GAMG.
> >
> [0] PCSetUp_GAMG(): level 0) N=120726, n data rows=3, n data cols=6,
> nnz/row (ave)=41, np=1
> [0] PCGAMGFilterGraph(): 99.904% nnz after filtering, with
> threshold 0., 13.7468 nnz ave. (N=40242)
> [0] PCGAMGCoarsen_AGG(): Square Graph on level 1 of 1 to square
> [0] PCGAMGProlongator_AGG(): New grid 1894 nodes
> [0] PCGAMGOptProlongator_AGG(): Smooth P0: max eigen=5.726852e+00
> min=1.330683e-01 PC=jacobi
> [0] PCSetUp_GAMG(): 1) N=11364, n data cols=6, nnz/row (ave)=196, 1
> active pes
> [0] PCGAMGFilterGraph(): 99.9839% nnz after filtering, with
> threshold 0., 32.7656 nnz ave. (N=1894)
> [0] PCGAMGProlongator_AGG(): New grid 155 nodes
> [0] PCGAMGOptProlongator_AGG(): Smooth P0: max eigen=1.09e+01
> min=1.832878e-04 PC=jacobi
> [0] PCSetUp_GAMG(): 2) N=930, n data cols=6, nnz/row (ave)=196, 1 active
> pes
> [0] PCGAMGFilterGraph(): 100.% nnz after filtering, with
> threshold 0., 32.7806 nnz ave. (N=155)
> [0] PCGAMGProlongator_AGG(): New grid 9 nodes
> [0] PCGAMGOptProlongator_AGG(): Smooth P0: max eigen=2.116373e+00
> min=6.337173e-03 PC=jacobi
> [0] PCSetUp_GAMG(): 3) N=54, n data cols=6, nnz/row (ave)=34, 1 active
> pes
> [0] PCGAMGFilterGraph(): 100.% nnz after filtering, with
> threshold 0., 5.7 nnz ave. (N=9)
> [0] PCGAMGProlongator_AGG(): New grid 2 nodes
> [0] PCGAMGOptProlongator_AGG(): Smooth P0: max eigen=1.984549e+00
> min=8.582767e-03 PC=jacobi
> [0] PCSetUp_GAMG(): 4) N=12, n data cols=6, nnz/row (ave)=12, 1 active
> pes
> [0] PCSetUp_GAMG(): 5 levels, grid complexity = 1.48586
> error: PETSc error 77-0 'Eigen estimator failed: DIVERGED_NANORINF at
> iteration 0'
> in /home/gtheler/libs/petsc-3.7.4/src/ksp/ksp/impls/cheby/cheby.c
> KSPSolve_Chebyshev:440
>
>
>
> >
> >
> >
> >
>
>
>
>
>

Re: [petsc-users] GAMG

[Jeremy Theler]

On Fri, 2016-10-28 at 09:46 -0400, Mark Adams wrote:
> Please run with -info and grep on GAMG.
> 
[0] PCSetUp_GAMG(): level 0) N=120726, n data rows=3, n data cols=6,
nnz/row (ave)=41, np=1
[0] PCGAMGFilterGraph(): 99.904% nnz after filtering, with
threshold 0., 13.7468 nnz ave. (N=40242)
[0] PCGAMGCoarsen_AGG(): Square Graph on level 1 of 1 to square
[0] PCGAMGProlongator_AGG(): New grid 1894 nodes
[0] PCGAMGOptProlongator_AGG(): Smooth P0: max eigen=5.726852e+00
min=1.330683e-01 PC=jacobi
[0] PCSetUp_GAMG(): 1) N=11364, n data cols=6, nnz/row (ave)=196, 1
active pes
[0] PCGAMGFilterGraph(): 99.9839% nnz after filtering, with
threshold 0., 32.7656 nnz ave. (N=1894)
[0] PCGAMGProlongator_AGG(): New grid 155 nodes
[0] PCGAMGOptProlongator_AGG(): Smooth P0: max eigen=1.09e+01
min=1.832878e-04 PC=jacobi
[0] PCSetUp_GAMG(): 2) N=930, n data cols=6, nnz/row (ave)=196, 1 active
pes
[0] PCGAMGFilterGraph(): 100.% nnz after filtering, with
threshold 0., 32.7806 nnz ave. (N=155)
[0] PCGAMGProlongator_AGG(): New grid 9 nodes
[0] PCGAMGOptProlongator_AGG(): Smooth P0: max eigen=2.116373e+00
min=6.337173e-03 PC=jacobi
[0] PCSetUp_GAMG(): 3) N=54, n data cols=6, nnz/row (ave)=34, 1 active
pes
[0] PCGAMGFilterGraph(): 100.% nnz after filtering, with
threshold 0., 5.7 nnz ave. (N=9)
[0] PCGAMGProlongator_AGG(): New grid 2 nodes
[0] PCGAMGOptProlongator_AGG(): Smooth P0: max eigen=1.984549e+00
min=8.582767e-03 PC=jacobi
[0] PCSetUp_GAMG(): 4) N=12, n data cols=6, nnz/row (ave)=12, 1 active
pes
[0] PCSetUp_GAMG(): 5 levels, grid complexity = 1.48586
error: PETSc error 77-0 'Eigen estimator failed: DIVERGED_NANORINF at
iteration 0'
in /home/gtheler/libs/petsc-3.7.4/src/ksp/ksp/impls/cheby/cheby.c
KSPSolve_Chebyshev:440



> 
> 
> 
> 

Re: [petsc-users] GAMG

[Jeremy Theler]

> 
> If I do not call PCSetCoordinates() the error goes away but
> convergence
> is slow.
> Is it possible that your coordinates lie on a 2D surface? All this
> does is make the 6 basis vectors
> for translations and rotations. You can just make these yourself and
> call MatSetNearNullSpace()
> and see what you get.
> 
No, they do not lie on a 2D surface :-/

Sorry but I did not get the point about the 6 basis vectors and
MatSetNearNullSpace().

--
jeremy

Re: [petsc-users] GAMG

[Jeremy Theler]

Hi Mark.

The matrix is solved well with lu/preonly.

If I do not call PCSetCoordinates() the error goes away but convergence
is slow.

I call PCSetCoordinates() this way (1 processor):

PetscMalloc1(dimensions * mesh->n_nodes, );
for (j = 0; j < mesh->n_nodes; j++) {
  for (d = 0; d < dimensions; d++) {
coords[j*dimensions + d] = mesh->node[j].x[d];
  }
}
PCSetCoordinates(pc, dimensions, dimensions * mesh->n_nodes,
coords);
PetscFree(coords);


Thanks
--
jeremy



On Fri, 2016-10-28 at 09:16 -0400, Mark Adams wrote:
> I think there is something wrong with your matrix. Use any solver and
> verify that you like the solution first.
> 
> On Fri, Oct 28, 2016 at 9:13 AM, Jeremy Theler <jer...@seamplex.com>
> wrote:
> Hi! I want to use PCGAMG as a preconditioner for a 3D linear
> elasticity
> problem (displacement-based FEM formulation) over an
> unstructured grid.
> I am not using DMPlex, I just build the stiffness matrix
> myself and pass
> it to PETSc.
> 
> I set MatSetBlockSize() to 3 and pass the node coordinates
> through
> PCSetCoordinates(). But using gamg and gmres I get:
> 
> PETSc error 77-0 'Eigen estimator failed: DIVERGED_NANORINF at
> iteration
> 0'
> in /home/gtheler/libs/petsc-3.7.4/src/ksp/ksp/impls/cheby/cheby.c
> KSPSolve_Chebyshev:440
> 
> Any suggestion? Another PC/KSP combination to try?
> 
> Thanks
> --
> jeremy
> 
> 
> 
> 

[petsc-users] GAMG

[Jeremy Theler]

Hi! I want to use PCGAMG as a preconditioner for a 3D linear elasticity
problem (displacement-based FEM formulation) over an unstructured grid.
I am not using DMPlex, I just build the stiffness matrix myself and pass
it to PETSc.

I set MatSetBlockSize() to 3 and pass the node coordinates through
PCSetCoordinates(). But using gamg and gmres I get:

PETSc error 77-0 'Eigen estimator failed: DIVERGED_NANORINF at iteration
0' in /home/gtheler/libs/petsc-3.7.4/src/ksp/ksp/impls/cheby/cheby.c
KSPSolve_Chebyshev:440

Any suggestion? Another PC/KSP combination to try?

Thanks
--
jeremy

Re: [petsc-users] Equivalent to MatGetColumnVector for rows?

[Jeremy Theler]

Thank you Barry. That makes a lot of sense.

--
jeremy

On Wed, 2016-10-19 at 14:54 -0500, Barry Smith wrote:
>I don't think you want to store the values in a Vector; the vector will be 
> as large as the entire right hand side but be almost all zeros.
> 
>If you want to remember "the part of the matrix that is zeroed out by 
> MatZeroRows()" you can use MatGetSubMatrix() and request just the zeroed rows 
> but all the columns. This matrix will be in parallel, on each process it will 
> just have the "zero rows" for that process. If you multiply this matrix by 
> the solution vector you will get a "short" vector that on each process 
> contains the "reaction" for each each of the "removed row" on that process.
> 
>    Easy to implement.
> 
>Barry
> 
> > On Oct 19, 2016, at 4:38 AM, Jeremy Theler <jer...@seamplex.com> wrote:
> > 
> > Hi all
> > 
> > Is there an equivalent to MatGetColumnVector() but for getting rows of a
> > matrix as a vector?
> > 
> > What I want to do is to compute the reactions of the nodes that belong
> > to a Dirichlet boundary condition in a FEM linear elastic problem. I set
> > these BCs with MatZeroRows() with a one in the diagonal and the desired
> > displacement in the RHS vector. But before calling MatZeroRows(), I want
> > to “remember” what the row looked like so after solving the problem, if
> > I multipliy this original row by the solution vector I get the reaction
> > corresponding to that row's DOF.
> > 
> > I have implemented something with MatGetRow() that seems to work but it
> > is some lame I am even embarrased of sharing with the list what I have
> > done.
> > 
> > Any suggestion is welcome.
> > 
> > Thanks
> > --
> > jeremy theler
> > www.seamplex.com
> > 
> > 
> 

[petsc-users] Equivalent to MatGetColumnVector for rows?

[Jeremy Theler]

Hi all

Is there an equivalent to MatGetColumnVector() but for getting rows of a
matrix as a vector?

What I want to do is to compute the reactions of the nodes that belong
to a Dirichlet boundary condition in a FEM linear elastic problem. I set
these BCs with MatZeroRows() with a one in the diagonal and the desired
displacement in the RHS vector. But before calling MatZeroRows(), I want
to “remember” what the row looked like so after solving the problem, if
I multipliy this original row by the solution vector I get the reaction
corresponding to that row's DOF.

I have implemented something with MatGetRow() that seems to work but it
is some lame I am even embarrased of sharing with the list what I have
done.

Any suggestion is welcome.

Thanks
--
jeremy theler
www.seamplex.com

Re: [petsc-users] Autoconf tests

[jeremy theler]

I once made a quick hack, maybe you can start your dirty work from here

https://bitbucket.org/wasora/wasora/src/5a88abbac1a846f2a6ed0b4e585f6f3c0fedf2e7/m4/petsc.m4?at=default=file-view-default

--
jeremy

On Tue, Oct 11, 2016 at 5:31 PM Barry Smith  wrote:

>
>You don't want to get the debug mode from PETSC_ARCH since there may
> not be a PETSC_ARCH (for PETSc --prefix installs) or because the user did
> not put the string in it.  You can check for the PETSC_USE_DEBUG symbol in
> the petscconf.h file by linking a C program against  and #if
> defined(PETSC_USE_DEBUG).
>
>
>Barry
>
> > On Oct 11, 2016, at 3:12 PM, CLAUS HELLMUTH WARNER HETZER <
> cl...@olemiss.edu> wrote:
> >
> > Hi everybody-
> >
> > Figured I’d ask this here before I go reinventing the wheel.
> >
> > I’m writing an autoconf installer (the standard Linux configure/make
> package) for an acoustic wave propagation modeling package that builds
> PETSc and SLEPc as part of the installation process.  I’d like to be able
> to test for instances of PETSc already being installed on the user’s
> machine and, if possible, whether they’re the debug version.  I know I can
> check for the existence of the PETSC_DIR environmental variable, and parse
> the PETSC_ARCH variable for “debug”, and I’ll do that as a first pass, but
> has anybody written any M4 tests that are more reliable than those (i.e.
> actually attempting to link to the libraries)?  I had one user who had the
> libraries installed in /usr/local/bin but didn’t have the environmental
> variables set in their profile, so the linker was confused and it took a
> while to figure out what was going weird with the install.
> >
> > If not, I guess I’ll be putting on my Autoconf gloves and getting my
> hands dirty.
> >
> > Thanks
> > -Claus Hetzer
> >
> > --
> > Claus Hetzer
> > Senior Research and Development Engineer
> > National Center for Physical Acoustics
> > The University of Mississippi
> > 145 Hill Drive
> > PO Box 1848
> > University, MS 38677
> > cl...@olemiss.edu
> >
> >
> >
> >
> >
>
>