Re: [petsc-users] Help with PETSc signal handling

[Jose E. Roman]

This SLEPc example uses PetscTime() to exit EPSSolve() after a given elapsed 
time.
http://slepc.upv.es/documentation/current/src/eps/examples/tutorials/ex29.c.html
Jose

> El 10 nov 2017, a las 5:09, Smith, Barry F.  escribió:
> 
> 
> 
>> On Nov 7, 2017, at 6:47 AM, Gard Spreemann  wrote:
>> 
>> On Tuesday 7 November 2017 07:35:36 CET Mark Adams wrote:
>>> PETSc's signal handler is for segvs, etc. I don't know the details but I
>>> don't think we care about external signals.
>> 
>> I see. I'll sketch what I'm trying to achieve in case someone can
>> think of another approach.
>> 
>> I have some long-running SLEPc eigenvalue computations, and I'd like
>> to have SLURM signal my program that its time limit is drawing
>> near. In that case, my problem would set a flag and before the next
>> iteration of the SLEPc eigenvalue solver it would give up and save the
>> eigenvalues it has so far managed to obtain.
>> 
>> The only workaround I can think of would be to let my program keep
>> track of its time limit on its own and check it at each iteration.
> 
>  Hmm, I would be more inclined to go with a model that avoided signals. 
> Signals are best avoided, unless absolutely necessary, plus you need to have 
> some external script/program to send in the signal requiring added complexity 
> to your work flow. 
> 
>   Why not have an input option of the time limit for the run and each 
> "iteration" or whatever (maybe in a post step function) check if you getting 
> close to the time-limit and then do what you need. If each "iteration" takes 
> a long time and you are afraid the program will end before it checks the next 
> time at the end of an iteration you could track record the time of each 
> iteration so your monitor routine knows how long a timestep takes and can 
> estimate if another timestep is possible within the current amount of time 
> you have available and if there is not enough time triggers the saving.
> 
>  Barry
> 
> 
> 
>> 
>> There is no intention for PETSc to have handling of user-defined
>> signals?
>> 
>> 
>> Thanks.
>> 
>> -- Gard
>> 
>> 
> 

Re: [petsc-users] Newton methods that converge all the time

[Smith, Barry F.]

  Henrik,

   Please describe in some detail how you are handling phase change. If you 
have if () tests of any sort in your FormFunction() or FormJacobian() this can 
kill Newton's method. If you are using "variable switching" this WILL kill 
Newtons' method. Are you monkeying with phase definitions in TSPostStep or with 
SNESLineSearchSetPostCheck(). This will also kill Newton's method.

  Barry


> On Nov 7, 2017, at 3:19 AM, Buesing, Henrik  
> wrote:
> 
> Dear all, 
>  
> I am solving a system of nonlinear, transient PDEs. I am using Newton’s 
> method in every time step to solve the nonlinear algebraic equations. Of 
> course, Newton’s method only converges if the initial guess is sufficiently 
> close to the solution. 
> 
> This is often not the case and Newton’s method diverges. Then, I reduce the 
> time step and try again. This can become prohibitively costly, if the time 
> steps get very small. I am thus looking for variants of Newton’s method, 
> which have a bigger convergence radius or ideally converge all the time. 
>  
> I tried out the pseudo-timestepping described in 
> http://www.mcs.anl.gov/petsc/petsc-current/src/ts/examples/tutorials/ex1f.F.html.
>  
> However, this does converge even worse. I am seeing breakdown when I have 
> phase changes (e.g. liquid to two-phase).
>  
> I was under the impression that pseudo-timestepping should converge better. 
> Thus, my question: 
> 
> Am I doing something wrong or is it possible that Newton’s method converges 
> and pseudo-timestepping does not?
> 
> Thank you for any insight on this. 
> 
> Henrik
> 
> 
> 
>  
> --
> Dipl.-Math. Henrik Büsing
> Institute for Applied Geophysics and Geothermal Energy
> E.ON Energy Research Center
> RWTH Aachen University
> --
> Mathieustr. 10|Tel +49 (0)241 80 49907
> 52074 Aachen, Germany |Fax +49 (0)241 80 49889
> --
> http://www.eonerc.rwth-aachen.de/GGE
> hbues...@eonerc.rwth-aachen.de
> --

Re: [petsc-users] Coloring of a finite volume unstructured mesh

[Smith, Barry F.]

  To use the PETSc coloring based Jacobian computer (which uses finite 
differences) you absolutely have to be able to provide the nonzero structure of 
the Jacobian. 

   Now once you provide the nonzero structure of the Jacobian the PETSc 
MatColoring routines can actually compute the coloring for you.

   So in other words you need not worry about the coloring at all, you just 
need to worry about providing the nonzero structure. Since you are using a 
finite volume method presumably all your coupling is between faces? In this 
case explicitly computing the nonzero structure of the Jacobian is probably 
pretty straightforward and you should just do it. 

  Barry


> On Nov 7, 2017, at 10:13 AM, Matthew Knepley  wrote:
> 
> On Tue, Nov 7, 2017 at 10:18 AM, SIERRA-AUSIN Javier 
>  wrote:
> Hi,
> 
> I would like to  ask you concerning the computation of the Jacobian matrix 
> via finite difference and coloring of the connectivity graph.
> I wonder whether it is possible or not to color the Jacobian matrix of a 
> given solver that evaluates the RHS with its associated connectivity in the 
> global indeces of my solver (not PETSc).
> As well, if it is possible to do this from an already partioned domain in 
> parallel.
> All of this is better explained in this post : 
> https://scicomp.stackexchange.com/questions/28209/linking-petsc-with-an-already-parallel-in-house-finite-volume-solver
>  
> 
> The simplest thing you can do is to use the finite-difference Jacobian action 
> (MatMFFD). This is setup automatically by SNES
> if you give a FormFunction pointer, but no FormJacobian routine. Just tell 
> the PETSc Vecs to use your ParMetis layout (by
> setting the local sizes), and it should run fine in SNES.
> 
> However, usually you need some kind of preconditioning. Thus you either have 
> to form the Jacobian or some approximation. If
> you cannot form an approximation, then you can use coloring. Once option is 
> to create a DMPlex with your mesh information.
> This can be done in parallel after you have already partitioned with ParMetis 
> (as long as you know the "overlap" of vertices, or
> adjacency of cells). Then the coloring can be done automatically using that 
> DM information. Otherwise, you will have to supply
> a coloring to the SNES.
> 
>   Thanks,
> 
>  Matt
>  
> Thanks in advance,
> 
> Javier.
>  
> 
> 
> 
> -- 
> 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] Help with PETSc signal handling

[Smith, Barry F.]

> On Nov 7, 2017, at 6:47 AM, Gard Spreemann  wrote:
> 
> On Tuesday 7 November 2017 07:35:36 CET Mark Adams wrote:
>> PETSc's signal handler is for segvs, etc. I don't know the details but I
>> don't think we care about external signals.
> 
> I see. I'll sketch what I'm trying to achieve in case someone can
> think of another approach.
> 
> I have some long-running SLEPc eigenvalue computations, and I'd like
> to have SLURM signal my program that its time limit is drawing
> near. In that case, my problem would set a flag and before the next
> iteration of the SLEPc eigenvalue solver it would give up and save the
> eigenvalues it has so far managed to obtain.
> 
> The only workaround I can think of would be to let my program keep
> track of its time limit on its own and check it at each iteration.

  Hmm, I would be more inclined to go with a model that avoided signals. 
Signals are best avoided, unless absolutely necessary, plus you need to have 
some external script/program to send in the signal requiring added complexity 
to your work flow. 

   Why not have an input option of the time limit for the run and each 
"iteration" or whatever (maybe in a post step function) check if you getting 
close to the time-limit and then do what you need. If each "iteration" takes a 
long time and you are afraid the program will end before it checks the next 
time at the end of an iteration you could track record the time of each 
iteration so your monitor routine knows how long a timestep takes and can 
estimate if another timestep is possible within the current amount of time you 
have available and if there is not enough time triggers the saving.

  Barry



> 
> There is no intention for PETSc to have handling of user-defined
> signals?
> 
> 
> Thanks.
> 
> -- Gard
> 
> 

Re: [petsc-users] Newton LS - better results on single processor

[Smith, Barry F.]

> On Nov 9, 2017, at 3:33 PM, zakaryah .  wrote:
> 
> Hi Stefano - when I referred to the iterations, I was trying to point out 
> that my method solves a series of nonlinear systems, with the solution to the 
> first problem being used to initialize the state vector for the second 
> problem, etc.  The reason I mentioned that was I thought perhaps I can expect 
> the residuals from single process solve to differ from the residuals from 
> multiprocess solve by a very small amount, say machine precision or the 
> tolerance of the KSP/SNES, that would be fine normally.  But, if there is a 
> possibility that those differences are somehow amplified by each of the 
> iterations (solution->initial state), that could explain what I see.
> 
> I agree that it is more likely that I have a bug in my code but I'm having 
> trouble finding it.

  Run a tiny problem on one and two processes with LU linear solver and the 
same mesh. So in the first case all values live on the first process and in the 
second the same first half live on one process and the second half on the 
second process. 

  Now track the values in the actual vectors and matrices. For example you can 
just put in VecView() and MatView() on all objects you pass into the solver and 
then put them in the SNESComputeFunction/Jacobian routines. Print both the 
vectors inputed to these routines and the vectors/matrices created in the 
routines. The output differences from the two runs should be small, determine 
when they significantly vary. This will tell you the likely location of the bug 
in your source code. (For example if certain values of the Jacobian differ)

  Good luck, I've done this plenty of times and if it is a "parallelization" 
bug this will help you find it much faster than guessing where the problem is 
and trying code inspect to find the bug.

  Barry

> 
> I ran a small problem with -pc_type redundant -redundant_pc_type lu, as you 
> suggested.  What I think is the relevant portion of the output is here (i.e. 
> there are small differences in the KSP residuals and SNES residuals):
> 
> -n 1, first "iteration" as described above:
> 
>   0 SNES Function norm 6.053565720454e-02 
>   
> 
> 0 KSP Residual norm 4.883115701982e-05
>   
> 
> 
> 0 KSP preconditioned resid norm 4.883115701982e-05 true resid norm 
> 6.053565720454e-02 ||r(i)||/||b|| 1.e+00  
>
> 
> 1 KSP Residual norm 8.173640409069e-20
>   
> 
> 
> 1 KSP preconditioned resid norm 8.173640409069e-20 true resid norm 
> 1.742143029296e-16 ||r(i)||/||b|| 2.877879104227e-15  
>
> 
>   Linear solve converged due to CONVERGED_RTOL iterations 1   
>   
> 
> 
>   Line search: Using full step: fnorm 6.053565720454e-02 gnorm 
> 2.735518570862e-07
>
> 
>   1 SNES Function norm 2.735518570862e-07 
>   
> 
> 
> 0 KSP Residual norm 1.298536630766e-10
>   
> 
> 
> 0 KSP preconditioned resid norm 1.298536630766e-10 true resid norm 
> 2.735518570862e-07 ||r(i)||/||b|| 1.e+00  
>
> 
> 1 KSP Residual norm 2.152782096751e-25
>   
> 
> 
> 1 KSP preconditioned resid norm 2.152782096751e-25 true resid norm 
> 4.75202641e-22 ||r(i)||/||b|| 1.738447420279e-15  
>
> 
>   Linear solve converged due to CONVERGED_RTOL iterations 1   
>   
> 
> 
>   Line search: Using full step: fnorm 2.735518570862e-07 gnorm 
> 

Re: [petsc-users] Newton LS - better results on single processor

[zakaryah .]

Hi Stefano - when I referred to the iterations, I was trying to point out
that my method solves a series of nonlinear systems, with the solution to
the first problem being used to initialize the state vector for the second
problem, etc.  The reason I mentioned that was I thought perhaps I can
expect the residuals from single process solve to differ from the residuals
from multiprocess solve by a very small amount, say machine precision or
the tolerance of the KSP/SNES, that would be fine normally.  But, if there
is a possibility that those differences are somehow amplified by each of
the iterations (solution->initial state), that could explain what I see.

I agree that it is more likely that I have a bug in my code but I'm having
trouble finding it.

I ran a small problem with -pc_type redundant -redundant_pc_type lu, as you
suggested.  What I think is the relevant portion of the output is here
(i.e. there are small differences in the KSP residuals and SNES residuals):

-n 1, first "iteration" as described above:

  0 SNES Function norm 6.053565720454e-02



0 KSP Residual norm 4.883115701982e-05



0 KSP preconditioned resid norm 4.883115701982e-05 true resid norm
6.053565720454e-02 ||r(i)||/||b|| 1.e+00


1 KSP Residual norm 8.173640409069e-20



1 KSP preconditioned resid norm 8.173640409069e-20 true resid norm
1.742143029296e-16 ||r(i)||/||b|| 2.877879104227e-15


  Linear solve converged due to CONVERGED_RTOL iterations 1



  Line search: Using full step: fnorm 6.053565720454e-02 gnorm
2.735518570862e-07


  1 SNES Function norm 2.735518570862e-07



0 KSP Residual norm 1.298536630766e-10



0 KSP preconditioned resid norm 1.298536630766e-10 true resid norm
2.735518570862e-07 ||r(i)||/||b|| 1.e+00


1 KSP Residual norm 2.152782096751e-25



1 KSP preconditioned resid norm 2.152782096751e-25 true resid norm
4.75202641e-22 ||r(i)||/||b|| 1.738447420279e-15


  Linear solve converged due to CONVERGED_RTOL iterations 1



  Line search: Using full step: fnorm 2.735518570862e-07 gnorm
1.917989238989e-17


  2 SNES Function norm 1.917989238989e-17



Nonlinear solve converged due to CONVERGED_FNORM_RELATIVE iterations 2


-n 2, first "iteration" as described above:

  0 SNES Function norm 6.053565720454e-02



0 KSP Residual norm 4.883115701982e-05



0 KSP preconditioned resid norm 4.883115701982e-05 true resid norm
6.053565720454e-02 ||r(i)||/||b|| 1.e+00


1 KSP Residual norm 1.007084240718e-19



1 KSP preconditioned resid norm 1.007084240718e-19 true resid norm
1.868472589717e-16 ||r(i)||/||b|| 3.086565300520e-15


  Linear solve converged due to CONVERGED_RTOL iterations 1



  Line search: Using full step: fnorm 6.053565720454e-02 gnorm
2.735518570379e-07


  1 SNES Function norm 2.735518570379e-07



0 KSP Residual norm 1.298536630342e-10



0 KSP preconditioned resid norm 1.298536630342e-10 true resid norm
2.735518570379e-07 ||r(i)||/||b|| 1.e+00


1 KSP Residual norm 1.885083482938e-25



1 KSP preconditioned resid norm 1.885083482938e-25 true resid norm
4.735707460766e-22 ||r(i)||/||b|| 1.731191852267e-15


  Linear solve converged due to CONVERGED_RTOL iterations 1



  Line search: Using full step: fnorm 2.735518570379e-07 gnorm
1.851472273258e-17


  2 SNES Function norm 1.851472273258e-17

-n 1, final "iteration":

  0 SNES Function norm 9.695669610792e+01



0 KSP Residual norm 7.898912593878e-03



0 KSP preconditioned resid norm 7.898912593878e-03 true resid norm
9.695669610792e+01 ||r(i)||/||b|| 1.e+00


1 KSP Residual norm 1.720960785852e-17



1 KSP preconditioned resid norm 1.720960785852e-17 true resid norm
1.23721391e-13 ||r(i)||/||b|| 1.275941911237e-15


  Linear solve converged due to CONVERGED_RTOL iterations 1



  Line search: Using full step: fnorm 9.695669610792e+01 gnorm
1.026572731653e-01


  1 SNES Function norm 1.026572731653e-01



0 KSP Residual norm 1.382450412926e-04



0 KSP preconditioned resid norm 1.382450412926e-04 true resid norm
1.026572731653e-01 ||r(i)||/||b|| 1.e+00


1 KSP Residual norm 5.018078565710e-20



1 KSP preconditioned resid norm 5.018078565710e-20 true resid norm
9.031463071676e-17 ||r(i)||/||b|| 8.797684560673e-16


  Linear solve converged due to CONVERGED_RTOL iterations 1



  Line search: Using full step: fnorm 1.026572731653e-01 gnorm
7.982937980399e-06


  2 SNES Function norm 7.982937980399e-06



0 KSP Residual norm 4.223898196692e-08



0 KSP preconditioned resid norm 4.223898196692e-08 true resid norm
7.982937980399e-06 ||r(i)||/||b|| 1.e+00


1 KSP Residual norm 1.038123933240e-22



1 KSP preconditioned resid norm 1.038123933240e-22 true resid norm
3.213931469966e-20 ||r(i)||/||b|| 4.026000800530e-15


  Linear solve converged due to CONVERGED_RTOL iterations 1



  Line search: Using full step: fnorm 

Re: [petsc-users] GAMG advice

[David Nolte]

Hi Mark,

thanks for clarifying.
When I wrote the initial question I had somehow overlooked the fact that
the GAMG standard smoother was Chebychev while ML uses SOR. All the
other comments concerning threshold etc were based on this mistake.

The following settings work quite well, of course LU is used on the
coarse level.

    -pc_type gamg
    -pc_gamg_type agg
    -pc_gamg_threshold 0.03
    -pc_gamg_square_graph 10        # no effect ?
    -pc_gamg_sym_graph
    -mg_levels_ksp_type richardson
    -mg_levels_pc_type sor

-pc_gamg_agg_nsmooths 0 does not seem to improve the convergence.

The ksp view now looks like this: (does this seem reasonable?)


KSP Object: 4 MPI processes
  type: fgmres
    GMRES: restart=30, using Classical (unmodified) Gram-Schmidt
Orthogonalization with no iterative refinement
    GMRES: happy breakdown tolerance 1e-30
  maximum iterations=1
  tolerances:  relative=1e-06, absolute=1e-50, divergence=1.
  right preconditioning
  using nonzero initial guess
  using UNPRECONDITIONED norm type for convergence test
PC Object: 4 MPI processes
  type: gamg
    MG: type is MULTIPLICATIVE, levels=5 cycles=v
  Cycles per PCApply=1
  Using Galerkin computed coarse grid matrices
  GAMG specific options
    Threshold for dropping small values from graph 0.03
    AGG specific options
  Symmetric graph true
  Coarse grid solver -- level ---
    KSP Object:    (mg_coarse_) 4 MPI processes
  type: preonly
  maximum iterations=1, initial guess is zero
  tolerances:  relative=1e-05, absolute=1e-50, divergence=1.
  left preconditioning
  using NONE norm type for convergence test
    PC Object:    (mg_coarse_) 4 MPI processes
  type: bjacobi
    block Jacobi: number of blocks = 4
    Local solve is same for all blocks, in the following KSP and PC
objects:
  KSP Object:  (mg_coarse_sub_)   1 MPI processes
    type: preonly
    maximum iterations=1, initial guess is zero
    tolerances:  relative=1e-05, absolute=1e-50, divergence=1.
    left preconditioning
    using NONE norm type for convergence test
  PC Object:  (mg_coarse_sub_)   1 MPI processes
    type: lu
  LU: out-of-place factorization
  tolerance for zero pivot 2.22045e-14
  using diagonal shift on blocks to prevent zero pivot [INBLOCKS]
  matrix ordering: nd
  factor fill ratio given 5., needed 1.
    Factored matrix follows:
  Mat Object:   1 MPI processes
    type: seqaij
    rows=38, cols=38
    package used to perform factorization: petsc
    total: nonzeros=1444, allocated nonzeros=1444
    total number of mallocs used during MatSetValues calls =0
  using I-node routines: found 8 nodes, limit used is 5
    linear system matrix = precond matrix:
    Mat Object: 1 MPI processes
  type: seqaij
  rows=38, cols=38
  total: nonzeros=1444, allocated nonzeros=1444
  total number of mallocs used during MatSetValues calls =0
    using I-node routines: found 8 nodes, limit used is 5
  linear system matrix = precond matrix:
  Mat Object:   4 MPI processes
    type: mpiaij
    rows=38, cols=38
    total: nonzeros=1444, allocated nonzeros=1444
    total number of mallocs used during MatSetValues calls =0
  using I-node (on process 0) routines: found 8 nodes, limit
used is 5
  Down solver (pre-smoother) on level 1 ---
    KSP Object:    (mg_levels_1_) 4 MPI processes
  type: richardson
    Richardson: damping factor=1.
  maximum iterations=2
  tolerances:  relative=1e-05, absolute=1e-50, divergence=1.
  left preconditioning
  using nonzero initial guess
  using NONE norm type for convergence test
    PC Object:    (mg_levels_1_) 4 MPI processes
  type: sor
    SOR: type = local_symmetric, iterations = 1, local iterations =
1, omega = 1.
  linear system matrix = precond matrix:
  Mat Object:   4 MPI processes
    type: mpiaij
    rows=168, cols=168
    total: nonzeros=19874, allocated nonzeros=19874
    total number of mallocs used during MatSetValues calls =0
  using I-node (on process 0) routines: found 17 nodes, limit
used is 5
  Up solver (post-smoother) same as down solver (pre-smoother)
  Down solver (pre-smoother) on level 2 ---
    KSP Object:    (mg_levels_2_) 4 MPI processes
  type: richardson
    Richardson: damping factor=1.
  maximum iterations=2
  tolerances:  relative=1e-05, absolute=1e-50, divergence=1.
  left preconditioning
  using nonzero initial guess
  using NONE norm type for convergence test
    PC Object:    (mg_levels_2_) 4 MPI processes
  type: 

Re: [petsc-users] unsorted local columns in 3.8?

[Hong]

Mark:

> OK, well, just go with the Linux machine for the regression test. I will
> keep trying to reproduce this on my Mac with an O build.
>

Valgrind error occurs on linux machines with g-build. I cannot merge this
branch to maint until the bug is fixed.
Hong

>
> On Wed, Nov 8, 2017 at 12:24 PM, Hong  wrote:
>
>> mpiexec -n 2 valgrind ./ex56 -cells 2,2,1 -max_conv_its 3
>> -petscspace_order 2 -snes_max_it 2 -ksp_max_it 100 -ksp_type cg -ksp_rtol
>> 1.e-11 -ksp_norm_type unpreconditioned -snes_rtol 1.e-10 -pc_type gamg
>> -pc_gamg_type agg -pc_gamg_agg_nsmooths 1 -pc_gamg_coarse_eq_limit 10
>> -pc_gamg_reuse_interpolation true -pc_gamg_square_graph 1
>> -pc_gamg_threshold 0.05 -pc_gamg_threshold_scale .0 -snes_converged_reason
>> -use_mat_nearnullspace true -mg_levels_ksp_max_it 1 -mg_levels_ksp_type
>> chebyshev -mg_levels_esteig_ksp_type cg -mg_levels_esteig_ksp_max_it 10
>> -mg_levels_ksp_chebyshev_esteig 0,0.05,0,1.05 -mg_levels_pc_type jacobi
>> -pc_gamg_mat_partitioning_type parmetis -mat_block_size 3 -run_type 1
>>
>> ==30976== Invalid read of size 16
>> ==30976==at 0x8550946: dswap_k_NEHALEM (in
>> /usr/lib/openblas-base/libblas.so.3)
>> ==30976==by 0x7C6797F: dswap_ (in /usr/lib/openblas-base/libblas
>> .so.3)
>> ==30976==by 0x75B33B2: dgetri_ (in /usr/lib/lapack/liblapack.so.3.0)
>> ==30976==by 0x5E3CA5C: PetscFESetUp_Basic (dtfe.c:4012)
>> ==30976==by 0x5E320C9: PetscFESetUp (dtfe.c:3274)
>> ==30976==by 0x5E5786F: PetscFECreateDefault (dtfe.c:6749)
>> ==30976==by 0x41056E: main (ex56.c:395)
>> ==30976==  Address 0xdc650d0 is 52,480 bytes inside a block of size
>> 52,488 alloc'd
>> ==30976==at 0x4C2D110: memalign (in /usr/lib/valgrind/vgpreload_me
>> mcheck-amd64-linux.so)
>> ==30976==by 0x51590F6: PetscMallocAlign (mal.c:39)
>> ==30976==by 0x5E3C169: PetscFESetUp_Basic (dtfe.c:3983)
>> ==30976==by 0x5E320C9: PetscFESetUp (dtfe.c:3274)
>> ==30976==by 0x5E5786F: PetscFECreateDefault (dtfe.c:6749)
>> ==30976==by 0x41056E: main (ex56.c:395)
>>
>> You can fix it on branch hzhang/fix-submat_samerowdist.
>>
>> Hong
>>
>>
>> On Wed, Nov 8, 2017 at 11:01 AM, Mark Adams  wrote:
>>
>>>
>>>
>>> On Wed, Nov 8, 2017 at 11:09 AM, Hong  wrote:
>>>
 Mark:

> Hong, is
> >   0-cells: 12 12 0 0
> >   1-cells: 20 20 0 0
> >   2-cells: 11 11 0 0
> >   3-cells: 2 2 0 0
>
> from the old version?
>
 In O-build on my macPro, I get the above. In g-build, I get

>
>   0-cells: 8 8 8 8
>   1-cells: 12 12 12 12
>   2-cells: 6 6 6 6
>   3-cells: 1 1 1 1
>
 I get this on linux machine.
 Do you know why?

>>>
>>> I can not reproduce your O output. I will look at it later.
>>>
>>> Valgrind is failing on me right now. I will look into it but can you
>>> valgrind it?
>>>
>>>

 Hong

>
> On Tue, Nov 7, 2017 at 10:13 PM, Hong  wrote:
>
>> Mark:
>> I removed option '-ex56_dm_view'.
>> Hong
>>
>> Humm, this looks a little odd, but it may be OK.  Is this this
>>> diffing with the old non-repartition data? (more below)
>>>
>>> On Tue, Nov 7, 2017 at 11:45 AM, Hong  wrote:
>>>
 Mark,
 The fix is merged to next branch for tests which show diff as

 *** Testing: testexamples_PARMETIS ***
 5c5
 <   1 SNES Function norm 1.983e-10
 ---
 >   1 SNES Function norm 1.990e-10
 10,13c10,13
 <   0-cells: 8 8 8 8
 <   1-cells: 12 12 12 12
 <   2-cells: 6 6 6 6
 <   3-cells: 1 1 1 1


>>> I assume this is the old.
>>>
>>>
 ---
 >   0-cells: 12 12 0 0
 >   1-cells: 20 20 0 0
 >   2-cells: 11 11 0 0
 >   3-cells: 2 2 0 0
 15,18c15,18


>>> and this is the new.
>>>
>>> This is funny because the processors are not fully populated. This
>>> can happen on coarse grids and indeed it should happen in a test with 
>>> good
>>> coverage.
>>>
>>> I assume these diffs are views from coarse grids? That is, in the
>>> raw output files do you see fully populated fine grids, with no diffs, 
>>> and
>>> then the diffs come on coarse grids.
>>>
>>> Repartitioning the coarse grids can change the coarsening, It is
>>> possible that repartitioning causes faster coarsening (it does a little)
>>> and this faster coarsening is tripping the aggregation switch, which 
>>> gives
>>> us empty processors.
>>>
>>> Am I understanding this correctly ...
>>>
>>> Thanks,
>>> Mark
>>>
>>>
 <   boundary: 1 strata with value/size (1 (23))
 <   Face Sets: 4 strata with value/size (1 (1), 2 (1), 4 (1), 6 (1))
 <   marker: 1 strata with value/size (1 (15))
 <   depth: 4 

Re: [petsc-users] DMPlexVecGetClosure for DM forest

[Matthew Knepley]

On Thu, Nov 9, 2017 at 1:20 PM, Yann Jobic  wrote:

> Hello,
> I'm trying to access to the values of a p4est forest.
> I know how to do that by converting my forest to a DMPlex, and then use
> DMPlexVecGetClosure over the converted DM.
> However, i want to assign a label and access to the values of the forest
> directly.
>

What do you mean by "directly". p4est only has topology. We use a Section
to map points to values, just like Plex.

 Thanks,

Matt


> Is it possible ?
> Thanks,
> Yann
>
>
> ---
> L'absence de virus dans ce courrier électronique a été vérifiée par le
> logiciel antivirus Avast.
> https://www.avast.com/antivirus
>
>


-- 
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] DMPlexVecGetClosure for DM forest

[Yann Jobic]

Hello,
I'm trying to access to the values of a p4est forest.
I know how to do that by converting my forest to a DMPlex, and then use
DMPlexVecGetClosure over the converted DM.
However, i want to assign a label and access to the values of the forest
directly.
Is it possible ?
Thanks,
Yann


---
L'absence de virus dans ce courrier électronique a été vérifiée par le logiciel 
antivirus Avast.
https://www.avast.com/antivirus

Re: [petsc-users] Newton LS - better results on single processor

[zakaryah .]

Thanks Stefano, I will try what you suggest.

​Matt - my DM is a composite between the redundant field (loading
coefficient, which is included in the Newton solve in Riks' method) and the
displacements, which are represented by a 3D DA with 3 dof.  I am using
finite difference.

Probably my problem comes from confusion over how the composite DM is
organized.  I am using FormFunction()​, and within that I call
DMCompositeGetLocalVectors(), DMCompositeScatter(), DMDAVecGetArray(), and
for the Jacobian, DMCompositeGetLocalISs() and MatGetLocalSubmatrix() to
split J into Jbb, Jbh, Jhb, and Jhh, where b is the loading coefficient,
and h is the displacements).  The values of each submatrix are set using
MatSetValuesLocal().

​I'm most suspicious of the part of the Jacobian routine where I calculate
the rows of Jhb, the columns of Jbh, and the corresponding values.  I take
the DA coordinates and ix,iy,iz, then calculate the row of Jhb as
iz-info->gzs)*info->gym
+ (iy-info->gys))*info->gxm + (ix-info->gxs))*info->dof+c), where info is
the DA local info and c is the degree of freedom.  The same calculation is
performed for the column of Jbh.  I suspect that the indexing of the DA
vector is not so simple, but I don't know for a fact that I'm doing this
incorrectly nor how to do this properly.

​Thanks for all the help!​


On Nov 9, 2017 8:44 AM, "Matthew Knepley"  wrote:

On Thu, Nov 9, 2017 at 12:14 AM, zakaryah .  wrote:

> Well the saga of my problem continues.  As I described previously in an
> epic thread, I'm using the SNES to solve problems involving an elastic
> material on a rectangular grid, subjected to external forces.  In any case,
> I'm occasionally getting poor convergence using Newton's method with line
> search.  In troubleshooting by visualizing the residual, I saw that in data
> sets which had good convergence, the residual was nevertheless
> significantly larger along the boundary between different processors.
> Likewise, in data sets with poor convergence, the residual became very
> large on the boundary between different processors.  The residual is not
> significantly larger on the physical boundary, i.e. the global boundary.
> When I run on a single process, convergence seems to be good on all data
> sets.
>
> Any clues to fix this?
>

It sounds like something is wrong with communication across domains:

 - If this is FEM, it sounds like you are not adding contributions from the
other domain to shared vertices/edges/faces

 - If this is FDM/FVM, maybe the ghosts are not updated

What DM are you using? Are you using the Local assembly functions
(FormFunctionLocal), or just FormFunction()?

  Thanks,

 Matt

-- 
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] DMPlex distribution and loops over cells/faces for finite volumes

[Matthew Knepley]

On Thu, Nov 9, 2017 at 3:08 AM, Matteo Semplice 
wrote:

> Hi.
> I am using a DMPLex to store a grid (at present 2d symplicial but will
> need to work more in general), but after distributing it, I am finding it
> hard to organize my loops over the cells/faces.
>
> First things first: the mesh distribution. I do
>
>   DMPlexSetAdjacencyUseCone(dm, PETSC_TRUE);
>   DMPlexSetAdjacencyUseClosure(dm, PETSC_FALSE);
>

This gives you adjacent cells, but not their boundary. Usually what you
want for FV.


>   DMPlexDistribute(dm, 1, NULL, );
>
> The 1 in DMPlexDistribute is correct to give me 1 layer of ghost cells,
> right?
>

Yes.


> Using 2 instead of 1 should add more internal ghosts, right? Actually,
> this errors out with petsc3.7.7 from debian/stable... If you deem it worth,
> I will test again with petsc3.8.
>

It should not error. If possible, send me the mesh please. I have some
tests where 2 is successful.


> Secondly, looping over cells. I do
>
>   DMPlexGetHeightStratum(dm, 0, , );
>   DMPlexGetHybridBounds(dm, , NULL, NULL, NULL);
>   for (int c = cStart; c < cEnd; ++c) {
>   if ((cPhysical>=0)&&(c>=cPhysical))
> {code for ghost cells outside physical domain boundary}
>   else if ( ??? )
> {code for ghost cells}
>

I am not sure why you want to draw this distinction, but you can check to
see if the cell is present in the pointSF. If it is, then its not owned by
the process.
Here is the kind of FV loop I do (for example in TS ex11):


https://bitbucket.org/petsc/petsc/src/f9fed41fd6c28d171f9c644f97b15591be9df7d6/src/snes/utils/dmplexsnes.c?at=master=file-view-default#dmplexsnes.c-1634


>   else
> {code for cells owned by this process}
>   }
>
>   What should I use as a test for ghost cells outside processor boundary?
> I have seen that (on a 2d symplicial mesh) reading the ghost label with
> DMGetLabelValue(dm, "ghost", c, ) gives the value -1 for inner cells
> and domain boundary ghosts but 2 for processor boundary ghosts. Is this the
> correct test? What is the meaning of ghost=2? In general, should I test for
> any value >=0?
>

Yes, see the link I gave.


> Lastly, since I will do finite volumes, I'd like to loop over the faces to
> compute fluxes.
> DMPlexGetHeightStratum(dm, 1, , )
> gives me the start/endpoints for faces, but they include also faces
> between ghost cells that are therefore external to the cells owned by the
> processor. I see that for faces on processor boundary, ghost=-1 on one
> process and ghost=2 on the other process, which I like. However ghost=-1
> also on faces between ghosts cells and therefore to compute fluxes I should
> test for the ghost value of the face and for the ghost value of the cells
> on both sides... Is there a better way to loop over faces?
>

See the link.

 Thanks,

   Matt


> Thanks in advance for any suggestion.
>
> Matteo
>
>


-- 
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] SNES without SNESSetJacobian, snes_fd or snes_mf

[Jed Brown]

Florian Kauer  writes:

> Hi,
> what is the SNES solver actually doing when you do not provide a 
> jacobian and not explicitly select either finite differencing 
> approximation or matrix-free Newton-Krylov method?
>
> I just noticed that I mistakenly did this, but a good solution is found 
> anyway (and fast). So what is actually happening? Simple fixed-point 
> iteration?

Presumably you are using a DM that can provide a coloring, in which case
the sparse Jacobian is assembled using finite differencing with coloring.

[petsc-users] multi level octree AMR

[Yann JOBIC]

Hi,

I succeed in having a basic AMR running on my FE advection/diffusion 
problem ( https://mycore.core-cloud.net/index.php/s/DmiiTwKUpV9z5qL)


I now want to have multiple levels in my octree AMR from p4est.

I tried a lot with IS tagged arrays, but i didn't succeed so far. When i 
get from VecTaggerAbsoluteSetBox an IS array of the cells, it seems that 
the DM cell numbering is changed after multiple time iterations (from 
TSSolve), even if the DM is not changed at all. I tried to identify the 
cells to refine/coarse with ISDifference and ISExpand, from two 
different time solutions of the same DM. That made me wondering if i'm 
using the correct tool (IS array).


Am i in the right direction ?

Thanks,

Regards,

Yann

Re: [petsc-users] unsorted local columns in 3.8?

[Mark Adams]

OK, well, just go with the Linux machine for the regression test. I will
keep trying to reproduce this on my Mac with an O build.

On Wed, Nov 8, 2017 at 12:24 PM, Hong  wrote:

> mpiexec -n 2 valgrind ./ex56 -cells 2,2,1 -max_conv_its 3
> -petscspace_order 2 -snes_max_it 2 -ksp_max_it 100 -ksp_type cg -ksp_rtol
> 1.e-11 -ksp_norm_type unpreconditioned -snes_rtol 1.e-10 -pc_type gamg
> -pc_gamg_type agg -pc_gamg_agg_nsmooths 1 -pc_gamg_coarse_eq_limit 10
> -pc_gamg_reuse_interpolation true -pc_gamg_square_graph 1
> -pc_gamg_threshold 0.05 -pc_gamg_threshold_scale .0 -snes_converged_reason
> -use_mat_nearnullspace true -mg_levels_ksp_max_it 1 -mg_levels_ksp_type
> chebyshev -mg_levels_esteig_ksp_type cg -mg_levels_esteig_ksp_max_it 10
> -mg_levels_ksp_chebyshev_esteig 0,0.05,0,1.05 -mg_levels_pc_type jacobi
> -pc_gamg_mat_partitioning_type parmetis -mat_block_size 3 -run_type 1
>
> ==30976== Invalid read of size 16
> ==30976==at 0x8550946: dswap_k_NEHALEM (in /usr/lib/openblas-base/
> libblas.so.3)
> ==30976==by 0x7C6797F: dswap_ (in /usr/lib/openblas-base/libblas.so.3)
> ==30976==by 0x75B33B2: dgetri_ (in /usr/lib/lapack/liblapack.so.3.0)
> ==30976==by 0x5E3CA5C: PetscFESetUp_Basic (dtfe.c:4012)
> ==30976==by 0x5E320C9: PetscFESetUp (dtfe.c:3274)
> ==30976==by 0x5E5786F: PetscFECreateDefault (dtfe.c:6749)
> ==30976==by 0x41056E: main (ex56.c:395)
> ==30976==  Address 0xdc650d0 is 52,480 bytes inside a block of size 52,488
> alloc'd
> ==30976==at 0x4C2D110: memalign (in /usr/lib/valgrind/vgpreload_
> memcheck-amd64-linux.so)
> ==30976==by 0x51590F6: PetscMallocAlign (mal.c:39)
> ==30976==by 0x5E3C169: PetscFESetUp_Basic (dtfe.c:3983)
> ==30976==by 0x5E320C9: PetscFESetUp (dtfe.c:3274)
> ==30976==by 0x5E5786F: PetscFECreateDefault (dtfe.c:6749)
> ==30976==by 0x41056E: main (ex56.c:395)
>
> You can fix it on branch hzhang/fix-submat_samerowdist.
>
> Hong
>
>
> On Wed, Nov 8, 2017 at 11:01 AM, Mark Adams  wrote:
>
>>
>>
>> On Wed, Nov 8, 2017 at 11:09 AM, Hong  wrote:
>>
>>> Mark:
>>>
 Hong, is
 >   0-cells: 12 12 0 0
 >   1-cells: 20 20 0 0
 >   2-cells: 11 11 0 0
 >   3-cells: 2 2 0 0

 from the old version?

>>> In O-build on my macPro, I get the above. In g-build, I get
>>>

   0-cells: 8 8 8 8
   1-cells: 12 12 12 12
   2-cells: 6 6 6 6
   3-cells: 1 1 1 1

>>> I get this on linux machine.
>>> Do you know why?
>>>
>>
>> I can not reproduce your O output. I will look at it later.
>>
>> Valgrind is failing on me right now. I will look into it but can you
>> valgrind it?
>>
>>
>>>
>>> Hong
>>>

 On Tue, Nov 7, 2017 at 10:13 PM, Hong  wrote:

> Mark:
> I removed option '-ex56_dm_view'.
> Hong
>
> Humm, this looks a little odd, but it may be OK.  Is this this diffing
>> with the old non-repartition data? (more below)
>>
>> On Tue, Nov 7, 2017 at 11:45 AM, Hong  wrote:
>>
>>> Mark,
>>> The fix is merged to next branch for tests which show diff as
>>>
>>> *** Testing: testexamples_PARMETIS ***
>>> 5c5
>>> <   1 SNES Function norm 1.983e-10
>>> ---
>>> >   1 SNES Function norm 1.990e-10
>>> 10,13c10,13
>>> <   0-cells: 8 8 8 8
>>> <   1-cells: 12 12 12 12
>>> <   2-cells: 6 6 6 6
>>> <   3-cells: 1 1 1 1
>>>
>>>
>> I assume this is the old.
>>
>>
>>> ---
>>> >   0-cells: 12 12 0 0
>>> >   1-cells: 20 20 0 0
>>> >   2-cells: 11 11 0 0
>>> >   3-cells: 2 2 0 0
>>> 15,18c15,18
>>>
>>>
>> and this is the new.
>>
>> This is funny because the processors are not fully populated. This
>> can happen on coarse grids and indeed it should happen in a test with 
>> good
>> coverage.
>>
>> I assume these diffs are views from coarse grids? That is, in the raw
>> output files do you see fully populated fine grids, with no diffs, and 
>> then
>> the diffs come on coarse grids.
>>
>> Repartitioning the coarse grids can change the coarsening, It is
>> possible that repartitioning causes faster coarsening (it does a little)
>> and this faster coarsening is tripping the aggregation switch, which 
>> gives
>> us empty processors.
>>
>> Am I understanding this correctly ...
>>
>> Thanks,
>> Mark
>>
>>
>>> <   boundary: 1 strata with value/size (1 (23))
>>> <   Face Sets: 4 strata with value/size (1 (1), 2 (1), 4 (1), 6 (1))
>>> <   marker: 1 strata with value/size (1 (15))
>>> <   depth: 4 strata with value/size (0 (8), 1 (12), 2 (6), 3 (1))
>>> ---
>>> >   boundary: 1 strata with value/size (1 (39))
>>> >   Face Sets: 5 strata with value/size (1 (2), 2 (2), 4 (2), 5 (1), 6 
>>> > (1))
>>> >   marker: 1 strata with value/size (1 (27))
>>> >   

[petsc-users] SNES without SNESSetJacobian, snes_fd or snes_mf

[Florian Kauer]

Hi,
what is the SNES solver actually doing when you do not provide a 
jacobian and not explicitly select either finite differencing 
approximation or matrix-free Newton-Krylov method?


I just noticed that I mistakenly did this, but a good solution is found 
anyway (and fast). So what is actually happening? Simple fixed-point 
iteration?


Greetings,
Florian

Re: [petsc-users] Newton LS - better results on single processor

[Matthew Knepley]

On Thu, Nov 9, 2017 at 12:14 AM, zakaryah .  wrote:

> Well the saga of my problem continues.  As I described previously in an
> epic thread, I'm using the SNES to solve problems involving an elastic
> material on a rectangular grid, subjected to external forces.  In any case,
> I'm occasionally getting poor convergence using Newton's method with line
> search.  In troubleshooting by visualizing the residual, I saw that in data
> sets which had good convergence, the residual was nevertheless
> significantly larger along the boundary between different processors.
> Likewise, in data sets with poor convergence, the residual became very
> large on the boundary between different processors.  The residual is not
> significantly larger on the physical boundary, i.e. the global boundary.
> When I run on a single process, convergence seems to be good on all data
> sets.
>
> Any clues to fix this?
>

It sounds like something is wrong with communication across domains:

 - If this is FEM, it sounds like you are not adding contributions from the
other domain to shared vertices/edges/faces

 - If this is FDM/FVM, maybe the ghosts are not updated

What DM are you using? Are you using the Local assembly functions
(FormFunctionLocal), or just FormFunction()?

  Thanks,

 Matt

-- 
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] linking PETSC with SSL

[Lukas van de Wiel]

Hi Satish,

I am working my way through linking our real code to PETSc 3.8.1. The
improvement for integration with Fortran is great.

use petscksp

makes a world of difference.

With compliments and thanks!
Lukas

On 11/8/17, Satish Balay  wrote:
> Glad it works! Thanks for the update.
>
> Satish
>
> On Wed, 8 Nov 2017, Lukas van de Wiel wrote:
>
>> Ah, never mind, It works in 3.8.1!
>> Had to change my makefile also, of course... :-\
>>
>> Thanks and have a great day!
>>
>> Lukas
>>
>>
>> On 11/8/17, Lukas van de Wiel  wrote:
>> > Hi Satish,
>> >
>> > thank you for your quick reply!
>> >
>> > I have installed 3.8.1, but the problem persists.
>> > I have attached configure.log and test.log.
>> >
>> > Does that reveal anything?
>> >
>> > Cheers and thanks again
>> >
>> > Lukas
>> >
>> > On 11/8/17, Satish Balay  wrote:
>> >> Hm --with-ssl=0 should work. Can you do a fresh/clean build and see if
>> >> the
>> >> problem persists?
>> >>
>> >> If so send configure.log make.log and test.log for this build.
>> >>
>> >> BTW: 3.8 has better support for fortran usage - so you might consider
>> >> upgrading all the way to it.
>> >>
>> >> https://www.mcs.anl.gov/petsc/documentation/changes/38.html
>> >> http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Sys/UsingFortran.html#UsingFortran
>> >>
>> >> Satish
>> >>
>> >> On Wed, 8 Nov 2017, Lukas van de Wiel wrote:
>> >>
>> >>> Good day,
>> >>>
>> >>> during an upgrade of PETSc (from 3.4.2 to 3.7.7) and some changes
>> >>> needed to use MUMPs I ran into a problem with SSL that I could not
>> >>> really figure out.
>> >>>
>> >>> 
>> >>> I have configured PETSc 3.7.7 using: (explicitly note the option
>> >>> --with-ssl=0  )
>> >>>
>> >>>
>> >>> ./configure  \
>> >>> COPTFLAGS='-O3 -march=native -mtune=native' \
>> >>> CXXOPTFLAGS='-O3 -march=native -mtune=native' \
>> >>> FOPTFLAGS='-O3 -march=native -mtune=native' \
>> >>> --with-debugging=0 \
>> >>> --with-x=0 \
>> >>> --with-ssl=0 \
>> >>> --with-shared-libraries=0 \
>> >>> --download-metis \
>> >>> --download-parmetis \
>> >>> --download-fblaslapack \
>> >>> --download-scalapack \
>> >>> --download-openmpi \
>> >>> --download-mumps \
>> >>> --download-hypre \
>> >>> --download-ptscotch
>> >>>
>> >>>
>> >>> 
>> >>> Next I have a tiny Fortran program to illustrate the problem:
>> >>>
>> >>>
>> >>> program sslhuh
>> >>> implicit none
>> >>> #include "petsc/finclude/petscsys.h"
>> >>> write(*,*) "Hello World"
>> >>> end program
>> >>>
>> >>>
>> >>>
>> >>> 
>> >>> I compile it using makefile:
>> >>>
>> >>>
>> >>> petscDir = /net/home/gtecton/sw_dev/petsc-3.7.7
>> >>>
>> >>> all: compile link
>> >>>
>> >>> compile: sslhuh.F
>> >>> $(petscDir)/linux-gnu-x86_64/bin/mpif90 \
>> >>> -c \
>> >>> -o sslhuh.o \
>> >>> -std=f2008 \
>> >>> -ffree-form \
>> >>> -I$(petscDir)/include \
>> >>> -I$(petscDir)/include/petsc/finclude \
>> >>> sslhuh.F
>> >>>
>> >>> link: sslhuh.o
>> >>> $(petscDir)/linux-gnu-x86_64/bin/mpif90 \
>> >>> -o sslhuh \
>> >>> sslhuh.o \
>> >>> $(petscDir)/linux-gnu-x86_64/lib/libpetsc.a \
>> >>> $(petscDir)/linux-gnu-x86_64/lib/libHYPRE.a \
>> >>> $(petscDir)/linux-gnu-x86_64/lib/libzmumps.a \
>> >>> $(petscDir)/linux-gnu-x86_64/lib/libsmumps.a \
>> >>> $(petscDir)/linux-gnu-x86_64/lib/libdmumps.a \
>> >>> $(petscDir)/linux-gnu-x86_64/lib/libcmumps.a \
>> >>> $(petscDir)/linux-gnu-x86_64/lib/libmumps_common.a \
>> >>> $(petscDir)/linux-gnu-x86_64/lib/libesmumps.a \
>> >>> $(petscDir)/linux-gnu-x86_64/lib/libparmetis.a \
>> >>> $(petscDir)/linux-gnu-x86_64/lib/libmetis.a \
>> >>> $(petscDir)/linux-gnu-x86_64/lib/libpord.a \
>> >>> $(petscDir)/linux-gnu-x86_64/lib/libscalapack.a \
>> >>> $(petscDir)/linux-gnu-x86_64/lib/libflapack.a \
>> >>> $(petscDir)/linux-gnu-x86_64/lib/libfblas.a \
>> >>> $(petscDir)/linux-gnu-x86_64/lib/libptscotch.a \
>> >>> $(petscDir)/linux-gnu-x86_64/lib/libptscotcherr.a \
>> >>> $(petscDir)/linux-gnu-x86_64/lib/libptscotcherrexit.a \
>> >>> $(petscDir)/linux-gnu-x86_64/lib/libptscotchparmetis.a \
>> >>> $(petscDir)/linux-gnu-x86_64/lib/libscotch.a \
>> >>> $(petscDir)/linux-gnu-x86_64/lib/libscotcherr.a \
>> >>> $(petscDir)/linux-gnu-x86_64/lib/libscotcherrexit.a \
>> >>> $(petscDir)/linux-gnu-x86_64/lib/libscotchmetis.a \
>> >>> -ldl
>> >>>
>> >>>
>> >>> 
>> >>> Compiling works fine.
>> >>> Linking gives me the error:
>> >>>
>> >>>
>> >>> /net/home/gtecton/sw_dev/petsc-3.7.7/linux-gnu-x86_64/lib/libpetsc.a(client.o):
>> >>> In function `PetscSSLInitializeContext':
>> >>> client.c:(.text+0x3b): undefined reference to `SSLv23_method'
>> >>> 

[petsc-users] DMPlex distribution and loops over cells/faces for finite volumes

[Matteo Semplice]

Hi.
I am using a DMPLex to store a grid (at present 2d symplicial but will 
need to work more in general), but after distributing it, I am finding 
it hard to organize my loops over the cells/faces.


First things first: the mesh distribution. I do

  DMPlexSetAdjacencyUseCone(dm, PETSC_TRUE);
  DMPlexSetAdjacencyUseClosure(dm, PETSC_FALSE);
  DMPlexDistribute(dm, 1, NULL, );

The 1 in DMPlexDistribute is correct to give me 1 layer of ghost cells, 
right?
Using 2 instead of 1 should add more internal ghosts, right? Actually, 
this errors out with petsc3.7.7 from debian/stable... If you deem it 
worth, I will test again with petsc3.8.


Secondly, looping over cells. I do

  DMPlexGetHeightStratum(dm, 0, , );
  DMPlexGetHybridBounds(dm, , NULL, NULL, NULL);
  for (int c = cStart; c < cEnd; ++c) {
  if ((cPhysical>=0)&&(c>=cPhysical))
    {code for ghost cells outside physical domain boundary}
  else if ( ??? )
    {code for ghost cells}
  else
    {code for cells owned by this process}
  }

  What should I use as a test for ghost cells outside processor 
boundary? I have seen that (on a 2d symplicial mesh) reading the ghost 
label with DMGetLabelValue(dm, "ghost", c, ) gives the value -1 
for inner cells and domain boundary ghosts but 2 for processor boundary 
ghosts. Is this the correct test? What is the meaning of ghost=2? In 
general, should I test for any value >=0?


Lastly, since I will do finite volumes, I'd like to loop over the faces 
to compute fluxes.

DMPlexGetHeightStratum(dm, 1, , )
gives me the start/endpoints for faces, but they include also faces 
between ghost cells that are therefore external to the cells owned by 
the processor. I see that for faces on processor boundary, ghost=-1 on 
one process and ghost=2 on the other process, which I like. However 
ghost=-1 also on faces between ghosts cells and therefore to compute 
fluxes I should test for the ghost value of the face and for the ghost 
value of the cells on both sides... Is there a better way to loop over 
faces?


Thanks in advance for any suggestion.

    Matteo