date:20160301

Re: [petsc-users] Error - Out of memory. This could be due to allocating too large an object or bleeding by not properly ...

2016-03-01 Thread Barry Smith


> On Mar 1, 2016, at 10:19 PM, TAY wee-beng  wrote:
> 
> 
> On 26/2/2016 9:21 PM, Barry Smith wrote:
>>> On Feb 26, 2016, at 1:14 AM, TAY wee-beng 
>>>  wrote:
>>> 
>>> 
>>> On 26/2/2016 1:56 AM, Barry Smith wrote:
>>> 
Run a much smaller problem for a few time steps, making sure you free 
 all the objects at the end, with the option -malloc_dump this will print 
 all the memory that was not freed and hopefully help you track down which 
 objects you forgot to free.
 
   Barry
 
>>> Hi,
>>> 
>>> I run a smaller problem and lots of things are shown in the log. How can I 
>>> know which exactly are not freed from the memory?
>>> 
>>Everything in in the log represents unfreed memory. You need to hunt 
>> through all the objects you create and make sure you destroy all of them.
>> 
>>   Barry
>> 
> Hi,
> 
> I have some questions.
> 
> [0]Total space allocated 2274656 bytes
> [ 0]16 bytes PetscStrallocpy() line 188 in 
> /home/wtay/Codes/petsc-3.6.3/src/sys/utils/str.c
> [ 0]624 bytes ISLocalToGlobalMappingCreate() line 270 in 
> /home/wtay/Codes/petsc-3.6.3/src/vec/is/utils/isltog.c
> [ 0]16 bytes VecScatterCreateCommon_PtoS() line 2655 in 
> /home/wtay/Codes/petsc-3.6.3/src/vec/vec/utils/vpscat.c
> [ 0]16 bytes VecScatterCreateCommon_PtoS() line 2654 in 
> /home/wtay/Codes/petsc-3.6.3/src/vec/vec/utils/vpscat.c
> [ 0]1440 bytes VecScatterCreate_PtoS() line 2463 in 
> /home/wtay/Codes/petsc-3.6.3/src/vec/vec/utils/vpscat.c
> [ 0]1440 bytes VecScatterCreate_PtoS() line 2462 in /home/wtay/Codes/
> 
> 1. What does the [0] means? I get from [0] to [23]?

  It is the MPI process reporting the memory usage
> 
> 2. I defined a variable globally:
> 
> DM  da_cu_types
> 
> Then I use at each time step:
> 
> call 
> DMDACreate3d(MPI_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,(IIB_I_end_domain(1)
>  - IIB_I_sta_domain(1) + 1),(IIB_I_end_domain(2) - IIB_I_sta_domain(2) + 1),&
> 
> (IIB_I_end_domain(3) - IIB_I_sta_domain(3) + 
> 1),PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,1,stencil_width_IIB,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,da_cu_types,ierr)
> 
> call 
> DMDAGetInfo(da_cu_types,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,num_procs_xyz_IIB(1),num_procs_xyz_IIB(2),num_procs_xyz_IIB(3),&
> 
> PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr)
> 
> call 
> DMDAGetCorners(da_cu_types,start_ijk_IIB(1),start_ijk_IIB(2),start_ijk_IIB(3),width_ijk_IIB(1),width_ijk_IIB(2),width_ijk_IIB(3),ierr)
> 
> call 
> DMDAGetGhostCorners(da_cu_types,start_ijk_ghost_IIB(1),start_ijk_ghost_IIB(2),start_ijk_ghost_IIB(3),width_ijk_ghost_IIB(1),width_ijk_ghost_IIB(2),width_ijk_ghost_IIB(3),ierr)
> 
> The purpose is just to get the starting and ending inidices for each cpu 
> partition. This is done every time step for a moving body case since the 
> IIB_I_sta_domain and IIB_I_end_domain changes
> 
> After getting all the info, must I call DMDestroy(da_cu_types,ierr)?

  Yes, otherwise you will get more and more DM taking up memory
> 
> Is it possible to update and use the new IIB_I_sta_domain and 
> IIB_I_end_domain without the need to create and destroy the DM? I thought 
> that may save some time since it's done at every time step.

  If the information you pass to the DMCreate changes then you need to create 
it again.

  Barry

> 
> Thanks
> 
>> 
>>> Is this info helpful? Or should I run in a single core?
>>> 
>>> Thanks
>>> 
> On Feb 25, 2016, at 12:33 AM, TAY wee-beng 
>  wrote:
> 
> Hi,
> 
> I ran the code and it hangs again. However, adding -malloc_test doesn't 
> seem to do any thing. The output (attached) is the same w/o it.
> 
> Wonder if there's anything else I can do.
> Thank you
> 
> Yours sincerely,
> 
> TAY wee-beng
> 
> On 24/2/2016 11:33 PM, Matthew Knepley wrote:
> 
>> On Wed, Feb 24, 2016 at 9:28 AM, TAY wee-beng 
>>  wrote:
>> 
>> On 24/2/2016 11:18 PM, Matthew Knepley wrote:
>> 
>>> On Wed, Feb 24, 2016 at 9:16 AM, TAY wee-beng 
>>>  wrote:
>>> 
>>> On 24/2/2016 9:12 PM, Matthew Knepley wrote:
>>> 
 On Wed, Feb 24, 2016 at 1:54 AM, TAY wee-beng 
  wrote:
 
 On 24/2/2016 10:28 AM, Matthew Knepley wrote:
 
> On Tue, Feb 23, 2016 at 7:50 PM, TAY wee-beng 
>  wrote:
> Hi,
> 
> I got this error (also attached, full) when running my code. It 
> happens after a few thousand time steps.
> 
> The strange thing is that for 2 different clusters, it stops at 2 
> different time steps.
> 
> I wonder if it's related to DM since this happens after I added DM 
> into my code.
>

Re: [petsc-users] Error - Out of memory. This could be due to allocating too large an object or bleeding by not properly ...

2016-03-01 Thread TAY wee-beng



On 26/2/2016 9:21 PM, Barry Smith wrote:

On Feb 26, 2016, at 1:14 AM, TAY wee-beng  wrote:


On 26/2/2016 1:56 AM, Barry Smith wrote:

Run a much smaller problem for a few time steps, making sure you free all 
the objects at the end, with the option -malloc_dump this will print all the 
memory that was not freed and hopefully help you track down which objects you 
forgot to free.

   Barry

Hi,

I run a smaller problem and lots of things are shown in the log. How can I know 
which exactly are not freed from the memory?

Everything in in the log represents unfreed memory. You need to hunt 
through all the objects you create and make sure you destroy all of them.

   Barry

Hi,

I have some questions.

[0]Total space allocated 2274656 bytes
[ 0]16 bytes PetscStrallocpy() line 188 in 
/home/wtay/Codes/petsc-3.6.3/src/sys/utils/str.c
[ 0]624 bytes ISLocalToGlobalMappingCreate() line 270 in 
/home/wtay/Codes/petsc-3.6.3/src/vec/is/utils/isltog.c
[ 0]16 bytes VecScatterCreateCommon_PtoS() line 2655 in 
/home/wtay/Codes/petsc-3.6.3/src/vec/vec/utils/vpscat.c
[ 0]16 bytes VecScatterCreateCommon_PtoS() line 2654 in 
/home/wtay/Codes/petsc-3.6.3/src/vec/vec/utils/vpscat.c
[ 0]1440 bytes VecScatterCreate_PtoS() line 2463 in 
/home/wtay/Codes/petsc-3.6.3/src/vec/vec/utils/vpscat.c

[ 0]1440 bytes VecScatterCreate_PtoS() line 2462 in /home/wtay/Codes/

1. What does the [0] means? I get from [0] to [23]?

2. I defined a variable globally:

DM  da_cu_types

Then I use at each time step:

/call 
DMDACreate3d(MPI_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,(IIB_I_end_domain(1) 
- IIB_I_sta_domain(1) + 1),(IIB_I_end_domain(2) - IIB_I_sta_domain(2) + 
1),&//

//
//(IIB_I_end_domain(3) - IIB_I_sta_domain(3) + 
1),PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,1,stencil_width_IIB,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,da_cu_types,ierr)//

//
//call 
DMDAGetInfo(da_cu_types,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,num_procs_xyz_IIB(1),num_procs_xyz_IIB(2),num_procs_xyz_IIB(3),&//

//
//PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr)//
//
//call 
DMDAGetCorners(da_cu_types,start_ijk_IIB(1),start_ijk_IIB(2),start_ijk_IIB(3),width_ijk_IIB(1),width_ijk_IIB(2),width_ijk_IIB(3),ierr)//

//
//call 
DMDAGetGhostCorners(da_cu_types,start_ijk_ghost_IIB(1),start_ijk_ghost_IIB(2),start_ijk_ghost_IIB(3),width_ijk_ghost_IIB(1),width_ijk_ghost_IIB(2),width_ijk_ghost_IIB(3),ierr)/


The purpose is just to get the starting and ending inidices for each cpu 
partition. This is done every time step for a moving body case since the 
/IIB_I_sta_domain/ and /IIB_I_end_domain/ changes


After getting all the info, must I call DMDestroy(da_cu_types,ierr)?

Is it possible to update and use the new /IIB_I_sta_domain/ and 
/IIB_I_end_domain /without the need to create and destroy the DM? I 
thought that may save some time since it's done at every time step.


Thanks




Is this info helpful? Or should I run in a single core?

Thanks

On Feb 25, 2016, at 12:33 AM, TAY wee-beng  wrote:

Hi,

I ran the code and it hangs again. However, adding -malloc_test doesn't seem to 
do any thing. The output (attached) is the same w/o it.

Wonder if there's anything else I can do.
Thank you

Yours sincerely,

TAY wee-beng

On 24/2/2016 11:33 PM, Matthew Knepley wrote:

On Wed, Feb 24, 2016 at 9:28 AM, TAY wee-beng  wrote:

On 24/2/2016 11:18 PM, Matthew Knepley wrote:

On Wed, Feb 24, 2016 at 9:16 AM, TAY wee-beng  wrote:

On 24/2/2016 9:12 PM, Matthew Knepley wrote:

On Wed, Feb 24, 2016 at 1:54 AM, TAY wee-beng  wrote:

On 24/2/2016 10:28 AM, Matthew Knepley wrote:

On Tue, Feb 23, 2016 at 7:50 PM, TAY wee-beng  wrote:
Hi,

I got this error (also attached, full) when running my code. It happens after a 
few thousand time steps.

The strange thing is that for 2 different clusters, it stops at 2 different 
time steps.

I wonder if it's related to DM since this happens after I added DM into my code.

In this case, how can I find out the error? I'm thinking valgrind may take very 
long and gives too many false errors.

It is very easy to find leaks. You just run a few steps with -malloc_dump and 
see what is left over.

Matt

Hi Matt,

Do you mean running my a.out with the -malloc_dump and stop after a few time 
steps?

What and how should I "see" then?

-malloc_dump outputs all unfreed memory to the screen after PetscFinalize(), so 
you should see the leak.
I guess it might be possible to keep creating things that you freed all at once 
at the end, but that is less likely.

Matt
  

Hi,

I got the output. I have zipped it since it's rather big. So it seems to be 
from DM routines but can you help me where the error is from?

Its really hard to tell by looking at it. What I do is remove things

Re: [petsc-users] Neumann BC with non-symmetric matrix

2016-03-01 Thread Mohammad Mirzadeh

On Tue, Mar 1, 2016 at 3:03 PM, Boyce Griffith 
wrote:

>
> On Mar 1, 2016, at 2:41 PM, Mohammad Mirzadeh  wrote:
>
>
>
> On Tue, Mar 1, 2016 at 2:07 PM, Boyce Griffith 
> wrote:
>
>>
>> On Mar 1, 2016, at 12:06 PM, Mohammad Mirzadeh 
>> wrote:
>>
>> Nice discussion.
>>
>>
>> On Tue, Mar 1, 2016 at 10:16 AM, Boyce Griffith 
>> wrote:
>>
>>>
>>> On Mar 1, 2016, at 9:59 AM, Mark Adams  wrote:
>>>
>>>
>>>
>>> On Mon, Feb 29, 2016 at 5:42 PM, Boyce Griffith 
>>> wrote:
>>>

 On Feb 29, 2016, at 5:36 PM, Mark Adams  wrote:


>> GAMG is use for AMR problems like this a lot in BISICLES.
>>
>
> Thanks for the reference. However, a quick look at their paper
> suggests they are using a finite volume discretization which should be
> symmetric and avoid all the shenanigans I'm going through!
>

 No, they are not symmetric.  FV is even worse than vertex centered
 methods.  The BCs and the C-F interfaces add non-symmetry.


 If you use a different discretization, it is possible to make the c-f
 interface discretization symmetric --- but symmetry appears to come at a
 cost of the reduction in the formal order of accuracy in the flux along the
 c-f interface. I can probably dig up some code that would make it easy to
 compare.

>>>
>>> I don't know.  Chombo/Boxlib have a stencil for C-F and do F-C with
>>> refluxing, which I do not linearize.  PETSc sums fluxes at faces directly,
>>> perhaps this IS symmetric? Toby might know.
>>>
>>>
>>> If you are talking about solving Poisson on a composite grid, then
>>> refluxing and summing up fluxes are probably the same procedure.
>>>
>>
>> I am not familiar with the terminology used here. What does the refluxing
>> mean?
>>
>>
>>>
>>> Users of these kinds of discretizations usually want to use the
>>> conservative divergence at coarse-fine interfaces, and so the main question
>>> is how to set up the viscous/diffusive flux stencil at coarse-fine
>>> interfaces (or, equivalently, the stencil for evaluating ghost cell values
>>> at coarse-fine interfaces). It is possible to make the overall
>>> discretization symmetric if you use a particular stencil for the flux
>>> computation. I think this paper (
>>> http://www.ams.org/journals/mcom/1991-56-194/S0025-5718-1991-1066831-5/S0025-5718-1991-1066831-5.pdf)
>>> is one place to look. (This stuff is related to "mimetic finite difference"
>>> discretizations of Poisson.) This coarse-fine interface discretization
>>> winds up being symmetric (although possibly only w.r.t. a weighted inner
>>> product --- I can't remember the details), but the fluxes are only
>>> first-order accurate at coarse-fine interfaces.
>>>
>>>
>> Right. I think if the discretization is conservative, i.e. discretizing
>> div of grad, and is compact, i.e. only involves neighboring cells sharing a
>> common face, then it is possible to construct symmetric discretization. An
>> example, that I have used before in other contexts, is described here:
>> http://physbam.stanford.edu/~fedkiw/papers/stanford2004-02.pdf
>>
>> An interesting observation is although the fluxes are only first order
>> accurate, the final solution to the linear system exhibits super
>> convergence, i.e. second-order accurate, even in L_inf. Similar behavior is
>> observed with non-conservative, node-based finite difference
>> discretizations.
>>
>>
>> I don't know about that --- check out Table 1 in the paper you cite,
>> which seems to indicate first-order convergence in all norms.
>>
>
> Sorry my bad. That was the original work which was later extended in
> doi:10.1016/j.compfluid.2005.01.006
>  to second order (c.f.
> section 3.3) by using flux weighting in the traverse direction.
>
>
> I don't follow the argument about why it is a bad thing for the fine
> fluxes to have different values than the overlying coarse flux, but this
> probably works out to be the same as the Ewing discretization in certain
> cases (although possibly only in 2D with a refinement ratio of 2).
>
>
In general, you can show that by weighting fluxes in the traverse direction
(a.k.a using the coarse-grid flux), you can eliminates the leading-order
truncation error (which turns out to be O(1)) in that direction. It seems,
to me, this is the reason you go from 1st to 2nd second-order accuracy.
Unfortunately this is not pointed out in the paper and while I cannot speak
for the authors, I think what motivated them to use the coarse flux was to
keep the normal velocity continuous across C-F interface after the
projection step.

It looks like that throughout the AMR literature people have encountered
the super-convergence over and over again and while most of the
"justifications" makes sense in an intuitive way, I don't think I

Re: [petsc-users] PCFactorSetUpMatSolverPackage with SNES

2016-03-01 Thread Barry Smith


> On Mar 1, 2016, at 6:43 PM, David Knezevic  wrote:
> 
> Hi Barry,
> 
> I have added an error checker so the code will produce a useful error 
> message instead of just crashing here in the code.
> 
> This sounds helpful. I assume this is in the dev branch?

  No it should work in 3.6.1 

> I'm using 3.6.1, but I gather from this list that 3.7 will be out soon, so 
> I'll switch to that once it's available.
> 
>  
> You can only call this routine AFTER the matrix has been provided to the 
> KSP object. Normally with SNES this won't happen until after the solve has 
> started (hence is not easy for you to put in your parameters) but adding the 
> code I suggest below might work
> 
> 
> OK, makes sense. Unfortunately the change below didn't help, I still get the 
> same segfault. But it's not a big deal, since I went with Matt's suggestion 
> which seems to be working well.
> 
> Thanks,
> David
> 
> 
>  
> > On Mar 1, 2016, at 12:52 PM, David Knezevic  
> > wrote:
> >
> > Based on KSP ex52, I use PCFactorSetUpMatSolverPackage in the process of 
> > setting various MUMPS ictnl options. This works fine for me when I'm 
> > solving linear problems.
> >
> > I then wanted to use PCFactorSetUpMatSolverPackage with the PC from a SNES 
> > object. I tried to do this with the following code (after calling 
> > SNESCreate, SNESSetFunction, and SNESSetJacobian):
> >
> > KSP snes_ksp;
> > SNESGetKSP(snes, _ksp);
> KSPSetOperators(ksp,mat,pmat);
>  /* mat and pmat (which may be the same) are the Jacobian arguments to 
> SNESSetJacobian(), the matrices must exist and types set but they don't have 
> to have the correct Jacobian entries at this point (since the nonlinear 
> solver has not started you cannot know the Jacobian entries yet.)
> > PC snes_pc;
> > KSPGetPC(snes_ksp, _pc);
> PCSetType(snes_pc,PCLU);
> > PCFactorSetMatSolverPackage(snes_pc, MATSOLVERMUMPS);
> > PCFactorSetUpMatSolverPackage(snes_pc);
> 
> Let me know if this does not work and where it goes wrong.
> >
> > However, I get a segfault on the call to PCFactorSetUpMatSolverPackage in 
> > this case. I was wondering what I need to do to make this work?
> >
> > Note that I want to set the MUMPS ictnl parameters via code rather than via 
> > the commandline since sometimes MUMPS fails (e.g. with error -9 due to a 
> > workspace size that is too small) and I need to automatically re-run the 
> > solve with different ictnl values when this happens.
> >
> > Thanks,
> > David
> >
> >
> 
>

Re: [petsc-users] PCFactorSetUpMatSolverPackage with SNES

2016-03-01 Thread David Knezevic

Hi Barry,

I have added an error checker so the code will produce a useful error
> message instead of just crashing here in the code.
>

This sounds helpful. I assume this is in the dev branch? I'm using 3.6.1,
but I gather from this list that 3.7 will be out soon, so I'll switch to
that once it's available.



> You can only call this routine AFTER the matrix has been provided to
> the KSP object. Normally with SNES this won't happen until after the solve
> has started (hence is not easy for you to put in your parameters) but
> adding the code I suggest below might work
>


OK, makes sense. Unfortunately the change below didn't help, I still get
the same segfault. But it's not a big deal, since I went with Matt's
suggestion which seems to be working well.

Thanks,
David




> > On Mar 1, 2016, at 12:52 PM, David Knezevic 
> wrote:
> >
> > Based on KSP ex52, I use PCFactorSetUpMatSolverPackage in the process of
> setting various MUMPS ictnl options. This works fine for me when I'm
> solving linear problems.
> >
> > I then wanted to use PCFactorSetUpMatSolverPackage with the PC from a
> SNES object. I tried to do this with the following code (after calling
> SNESCreate, SNESSetFunction, and SNESSetJacobian):
> >
> > KSP snes_ksp;
> > SNESGetKSP(snes, _ksp);
> KSPSetOperators(ksp,mat,pmat);
>  /* mat and pmat (which may be the same) are the Jacobian arguments to
> SNESSetJacobian(), the matrices must exist and types set but they don't
> have to have the correct Jacobian entries at this point (since the
> nonlinear solver has not started you cannot know the Jacobian entries yet.)
> > PC snes_pc;
> > KSPGetPC(snes_ksp, _pc);
> PCSetType(snes_pc,PCLU);
> > PCFactorSetMatSolverPackage(snes_pc, MATSOLVERMUMPS);
> > PCFactorSetUpMatSolverPackage(snes_pc);
>
> Let me know if this does not work and where it goes wrong.
> >
> > However, I get a segfault on the call to PCFactorSetUpMatSolverPackage
> in this case. I was wondering what I need to do to make this work?
> >
> > Note that I want to set the MUMPS ictnl parameters via code rather than
> via the commandline since sometimes MUMPS fails (e.g. with error -9 due to
> a workspace size that is too small) and I need to automatically re-run the
> solve with different ictnl values when this happens.
> >
> > Thanks,
> > David
> >
> >
>
>

Re: [petsc-users] PCFactorSetUpMatSolverPackage with SNES

2016-03-01 Thread Barry Smith


   David,

I have added an error checker so the code will produce a useful error 
message instead of just crashing here in the code.

You can only call this routine AFTER the matrix has been provided to the 
KSP object. Normally with SNES this won't happen until after the solve has 
started (hence is not easy for you to put in your parameters) but adding the 
code I suggest below might work

> On Mar 1, 2016, at 12:52 PM, David Knezevic  
> wrote:
> 
> Based on KSP ex52, I use PCFactorSetUpMatSolverPackage in the process of 
> setting various MUMPS ictnl options. This works fine for me when I'm solving 
> linear problems.
> 
> I then wanted to use PCFactorSetUpMatSolverPackage with the PC from a SNES 
> object. I tried to do this with the following code (after calling SNESCreate, 
> SNESSetFunction, and SNESSetJacobian):
> 
> KSP snes_ksp;
> SNESGetKSP(snes, _ksp);
KSPSetOperators(ksp,mat,pmat);   
 /* mat and pmat (which may be the same) are the Jacobian arguments to 
SNESSetJacobian(), the matrices must exist and types set but they don't have to 
have the correct Jacobian entries at this point (since the nonlinear solver has 
not started you cannot know the Jacobian entries yet.)
> PC snes_pc;
> KSPGetPC(snes_ksp, _pc);
PCSetType(snes_pc,PCLU);
> PCFactorSetMatSolverPackage(snes_pc, MATSOLVERMUMPS);
> PCFactorSetUpMatSolverPackage(snes_pc);

Let me know if this does not work and where it goes wrong.
> 
> However, I get a segfault on the call to PCFactorSetUpMatSolverPackage in 
> this case. I was wondering what I need to do to make this work?
> 
> Note that I want to set the MUMPS ictnl parameters via code rather than via 
> the commandline since sometimes MUMPS fails (e.g. with error -9 due to a 
> workspace size that is too small) and I need to automatically re-run the 
> solve with different ictnl values when this happens.
> 
> Thanks,
> David
> 
>

Re: [petsc-users] PCFactorSetUpMatSolverPackage with SNES

2016-03-01 Thread David Knezevic

On Tue, Mar 1, 2016 at 1:56 PM, Matthew Knepley  wrote:

> On Tue, Mar 1, 2016 at 12:52 PM, David Knezevic <
> david.kneze...@akselos.com> wrote:
>
>> Based on KSP ex52, I use PCFactorSetUpMatSolverPackage in the process of
>> setting various MUMPS ictnl options. This works fine for me when I'm
>> solving linear problems.
>>
>> I then wanted to use PCFactorSetUpMatSolverPackage with the PC from a
>> SNES object. I tried to do this with the following code (after calling
>> SNESCreate, SNESSetFunction, and SNESSetJacobian):
>>
>> KSP snes_ksp;
>> SNESGetKSP(snes, _ksp);
>> PC snes_pc;
>> KSPGetPC(snes_ksp, _pc);
>> PCFactorSetMatSolverPackage(snes_pc, MATSOLVERMUMPS);
>> PCFactorSetUpMatSolverPackage(snes_pc);
>>
>> However, I get a segfault on the call to PCFactorSetUpMatSolverPackage in
>> this case. I was wondering what I need to do to make this work?
>>
>> Note that I want to set the MUMPS ictnl parameters via code rather than
>> via the commandline since sometimes MUMPS fails (e.g. with error -9 due to
>> a workspace size that is too small) and I need to automatically re-run the
>> solve with different ictnl values when this happens.
>>
>
> That is a good reason. However, I would organize this differently. I would
> still set the type from the command line.
> Then later in your code, after SNES is setup correctly, you get out the PC
> and reset the icntl values if you have a
> failure. Its very difficult to get the setup logic correct, and its one of
> the most error prone parts of PETSc. RIght now,
> the most reliable way to do things is to have all the information
> available up front in options.
>
> Someday, I want to write a small DFA piece in PETSc that can encode all
> this logic so that simple errors people
> make can be diagnosed early with nice error messages.
>
>   Thanks,
>
>  Matt
>

OK, thanks for the info, I'll go with that approach.

David

Re: [petsc-users] Neumann BC with non-symmetric matrix

2016-03-01 Thread Boyce Griffith


> On Mar 1, 2016, at 2:41 PM, Mohammad Mirzadeh  wrote:
> 
> 
> 
> On Tue, Mar 1, 2016 at 2:07 PM, Boyce Griffith  > wrote:
> 
>> On Mar 1, 2016, at 12:06 PM, Mohammad Mirzadeh > > wrote:
>> 
>> Nice discussion.
>> 
>> 
>> On Tue, Mar 1, 2016 at 10:16 AM, Boyce Griffith > > wrote:
>> 
>>> On Mar 1, 2016, at 9:59 AM, Mark Adams >> > wrote:
>>> 
>>> 
>>> 
>>> On Mon, Feb 29, 2016 at 5:42 PM, Boyce Griffith >> > wrote:
>>> 
 On Feb 29, 2016, at 5:36 PM, Mark Adams > wrote:
 
 
 GAMG is use for AMR problems like this a lot in BISICLES.
 
 Thanks for the reference. However, a quick look at their paper suggests 
 they are using a finite volume discretization which should be symmetric 
 and avoid all the shenanigans I'm going through! 
 
 No, they are not symmetric.  FV is even worse than vertex centered 
 methods.  The BCs and the C-F interfaces add non-symmetry.
>>> 
>>> 
>>> If you use a different discretization, it is possible to make the c-f 
>>> interface discretization symmetric --- but symmetry appears to come at a 
>>> cost of the reduction in the formal order of accuracy in the flux along the 
>>> c-f interface. I can probably dig up some code that would make it easy to 
>>> compare.
>>> 
>>> I don't know.  Chombo/Boxlib have a stencil for C-F and do F-C with 
>>> refluxing, which I do not linearize.  PETSc sums fluxes at faces directly, 
>>> perhaps this IS symmetric? Toby might know.
>> 
>> If you are talking about solving Poisson on a composite grid, then refluxing 
>> and summing up fluxes are probably the same procedure.
>> 
>> I am not familiar with the terminology used here. What does the refluxing 
>> mean?
>>  
>> 
>> Users of these kinds of discretizations usually want to use the conservative 
>> divergence at coarse-fine interfaces, and so the main question is how to set 
>> up the viscous/diffusive flux stencil at coarse-fine interfaces (or, 
>> equivalently, the stencil for evaluating ghost cell values at coarse-fine 
>> interfaces). It is possible to make the overall discretization symmetric if 
>> you use a particular stencil for the flux computation. I think this paper 
>> (http://www.ams.org/journals/mcom/1991-56-194/S0025-5718-1991-1066831-5/S0025-5718-1991-1066831-5.pdf
>>  
>> )
>>  is one place to look. (This stuff is related to "mimetic finite difference" 
>> discretizations of Poisson.) This coarse-fine interface discretization winds 
>> up being symmetric (although possibly only w.r.t. a weighted inner product 
>> --- I can't remember the details), but the fluxes are only first-order 
>> accurate at coarse-fine interfaces.
>> 
>> 
>> Right. I think if the discretization is conservative, i.e. discretizing div 
>> of grad, and is compact, i.e. only involves neighboring cells sharing a 
>> common face, then it is possible to construct symmetric discretization. An 
>> example, that I have used before in other contexts, is described here: 
>> http://physbam.stanford.edu/~fedkiw/papers/stanford2004-02.pdf 
>> 
>> 
>> An interesting observation is although the fluxes are only first order 
>> accurate, the final solution to the linear system exhibits super 
>> convergence, i.e. second-order accurate, even in L_inf. Similar behavior is 
>> observed with non-conservative, node-based finite difference 
>> discretizations. 
> 
> I don't know about that --- check out Table 1 in the paper you cite, which 
> seems to indicate first-order convergence in all norms.
> 
> Sorry my bad. That was the original work which was later extended in 
> doi:10.1016/j.compfluid.2005.01.006 
>  to second order (c.f. 
> section 3.3) by using flux weighting in the traverse direction.

I don't follow the argument about why it is a bad thing for the fine fluxes to 
have different values than the overlying coarse flux, but this probably works 
out to be the same as the Ewing discretization in certain cases (although 
possibly only in 2D with a refinement ratio of 2).

> A (discrete) Green's functions argument explains why one gets higher-order 
> convergence despite localized reductions in accuracy along the coarse-fine 
> interface --- it has to do with the fact that errors from individual grid 
> locations do not have that large of an effect on the solution, and these c-f 
> interface errors are concentrated along on a lower dimensional surface in the 
> domain.
> 
> This intuitively makes sense and in fact when you plot the error, you do

Re: [petsc-users] Neumann BC with non-symmetric matrix

2016-03-01 Thread Matthew Knepley

On Tue, Mar 1, 2016 at 1:48 PM, Mohammad Mirzadeh 
wrote:

> On Tue, Mar 1, 2016 at 1:15 PM, Jed Brown  wrote:
>
>> Mohammad Mirzadeh  writes:
>>
>> > I am not familiar with the terminology used here. What does the
>> refluxing
>> > mean?
>>
>> The Chombo/BoxLib family of methods evaluate fluxes between coarse grid
>> cells overlaying refined grids, then later visit the fine grids and
>> reevaluate those fluxes.  The correction needs to be propagated back to
>> the adjoining coarse grid cell to maintain conservation.  It's an
>> implementation detail that they call refluxing.
>>
>
> Thanks for clarification.
>
>
>>
>> > Right. I think if the discretization is conservative, i.e. discretizing
>> div
>> > of grad, and is compact, i.e. only involves neighboring cells sharing a
>> > common face, then it is possible to construct symmetric discretization.
>> An
>> > example, that I have used before in other contexts, is described here:
>> > http://physbam.stanford.edu/~fedkiw/papers/stanford2004-02.pdf
>>
>> It's unfortunate that this paper repeats some unfounded multigrid
>> slander and then basically claims to have uniform convergence using
>> incomplete Cholesky with CG.  In reality, incomplete Cholesky is
>> asymptotically no better than Jacobi.
>>
>> > An interesting observation is although the fluxes are only first order
>> > accurate, the final solution to the linear system exhibits super
>> > convergence, i.e. second-order accurate, even in L_inf.
>>
>> Perhaps for aligned coefficients; definitely not for unaligned
>> coefficients.
>>
>
> Could you elaborate what you mean by aligned/unaligned coefficients? Do
> you mean anisotropic diffusion coefficient?
>

Jed (I think) means coefficients where the variation is aligned to the
grid. For example, where coefficient jumps or large variation
happens on cell boundaries.

   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

Re: [petsc-users] Neumann BC with non-symmetric matrix

2016-03-01 Thread Mohammad Mirzadeh

On Tue, Mar 1, 2016 at 1:15 PM, Jed Brown  wrote:

> Mohammad Mirzadeh  writes:
>
> > I am not familiar with the terminology used here. What does the refluxing
> > mean?
>
> The Chombo/BoxLib family of methods evaluate fluxes between coarse grid
> cells overlaying refined grids, then later visit the fine grids and
> reevaluate those fluxes.  The correction needs to be propagated back to
> the adjoining coarse grid cell to maintain conservation.  It's an
> implementation detail that they call refluxing.
>

Thanks for clarification.


>
> > Right. I think if the discretization is conservative, i.e. discretizing
> div
> > of grad, and is compact, i.e. only involves neighboring cells sharing a
> > common face, then it is possible to construct symmetric discretization.
> An
> > example, that I have used before in other contexts, is described here:
> > http://physbam.stanford.edu/~fedkiw/papers/stanford2004-02.pdf
>
> It's unfortunate that this paper repeats some unfounded multigrid
> slander and then basically claims to have uniform convergence using
> incomplete Cholesky with CG.  In reality, incomplete Cholesky is
> asymptotically no better than Jacobi.
>
> > An interesting observation is although the fluxes are only first order
> > accurate, the final solution to the linear system exhibits super
> > convergence, i.e. second-order accurate, even in L_inf.
>
> Perhaps for aligned coefficients; definitely not for unaligned
> coefficients.
>

Could you elaborate what you mean by aligned/unaligned coefficients? Do you
mean anisotropic diffusion coefficient?

Re: [petsc-users] Neumann BC with non-symmetric matrix

2016-03-01 Thread Mohammad Mirzadeh

On Tue, Mar 1, 2016 at 2:07 PM, Boyce Griffith 
wrote:

>
> On Mar 1, 2016, at 12:06 PM, Mohammad Mirzadeh  wrote:
>
> Nice discussion.
>
>
> On Tue, Mar 1, 2016 at 10:16 AM, Boyce Griffith 
> wrote:
>
>>
>> On Mar 1, 2016, at 9:59 AM, Mark Adams  wrote:
>>
>>
>>
>> On Mon, Feb 29, 2016 at 5:42 PM, Boyce Griffith 
>> wrote:
>>
>>>
>>> On Feb 29, 2016, at 5:36 PM, Mark Adams  wrote:
>>>
>>>
> GAMG is use for AMR problems like this a lot in BISICLES.
>

 Thanks for the reference. However, a quick look at their paper suggests
 they are using a finite volume discretization which should be symmetric and
 avoid all the shenanigans I'm going through!

>>>
>>> No, they are not symmetric.  FV is even worse than vertex centered
>>> methods.  The BCs and the C-F interfaces add non-symmetry.
>>>
>>>
>>> If you use a different discretization, it is possible to make the c-f
>>> interface discretization symmetric --- but symmetry appears to come at a
>>> cost of the reduction in the formal order of accuracy in the flux along the
>>> c-f interface. I can probably dig up some code that would make it easy to
>>> compare.
>>>
>>
>> I don't know.  Chombo/Boxlib have a stencil for C-F and do F-C with
>> refluxing, which I do not linearize.  PETSc sums fluxes at faces directly,
>> perhaps this IS symmetric? Toby might know.
>>
>>
>> If you are talking about solving Poisson on a composite grid, then
>> refluxing and summing up fluxes are probably the same procedure.
>>
>
> I am not familiar with the terminology used here. What does the refluxing
> mean?
>
>
>>
>> Users of these kinds of discretizations usually want to use the
>> conservative divergence at coarse-fine interfaces, and so the main question
>> is how to set up the viscous/diffusive flux stencil at coarse-fine
>> interfaces (or, equivalently, the stencil for evaluating ghost cell values
>> at coarse-fine interfaces). It is possible to make the overall
>> discretization symmetric if you use a particular stencil for the flux
>> computation. I think this paper (
>> http://www.ams.org/journals/mcom/1991-56-194/S0025-5718-1991-1066831-5/S0025-5718-1991-1066831-5.pdf)
>> is one place to look. (This stuff is related to "mimetic finite difference"
>> discretizations of Poisson.) This coarse-fine interface discretization
>> winds up being symmetric (although possibly only w.r.t. a weighted inner
>> product --- I can't remember the details), but the fluxes are only
>> first-order accurate at coarse-fine interfaces.
>>
>>
> Right. I think if the discretization is conservative, i.e. discretizing
> div of grad, and is compact, i.e. only involves neighboring cells sharing a
> common face, then it is possible to construct symmetric discretization. An
> example, that I have used before in other contexts, is described here:
> http://physbam.stanford.edu/~fedkiw/papers/stanford2004-02.pdf
>
> An interesting observation is although the fluxes are only first order
> accurate, the final solution to the linear system exhibits super
> convergence, i.e. second-order accurate, even in L_inf. Similar behavior is
> observed with non-conservative, node-based finite difference
> discretizations.
>
>
> I don't know about that --- check out Table 1 in the paper you cite, which
> seems to indicate first-order convergence in all norms.
>

Sorry my bad. That was the original work which was later extended in
doi:10.1016/j.compfluid.2005.01.006
 to second order (c.f.
section 3.3) by using flux weighting in the traverse direction.


>
> The symmetric discretization in the Ewing paper is only slightly more
> complicated, but will give full 2nd-order accuracy in L-1 (and maybe also
> L-2 and L-infinity). One way to think about it is that you are using simple
> linear interpolation at coarse-fine interfaces (3-point interpolation in
> 2D, 4-point interpolation in 3D) using a stencil that is symmetric with
> respect to the center of the coarse grid cell.
>
>
I'll look into that paper. One can never get enough of ideas for C-F
treatments in AMR applications :).


> A (discrete) Green's functions argument explains why one gets higher-order
> convergence despite localized reductions in accuracy along the coarse-fine
> interface --- it has to do with the fact that errors from individual grid
> locations do not have that large of an effect on the solution, and these
> c-f interface errors are concentrated along on a lower dimensional surface
> in the domain.
>

This intuitively makes sense and in fact when you plot the error, you do
see spikes at the C-F interfaces. Do you know of a resource that does a
rigorous analysis of the C-F treatment on the solution error?


>
> -- Boyce
>

Re: [petsc-users] Neumann BC with non-symmetric matrix

2016-03-01 Thread Griffith, Boyce Eugene


On Mar 1, 2016, at 10:56 AM, Jed Brown 
> wrote:

Boyce Griffith > writes:
Jed, can you also do this for Stokes?  It seems like something like
RT0 is the right place to start.

See, for example, Arnold, Falk, and Winther's 2007 paper on mixed FEM
for elasticity with weakly imposed symmetry.  It's the usual H(div)
methodology and should apply equally well to Stokes.  I'm not aware of
any analysis or results of choosing quadrature to eliminate flux terms
in these discretizations.

Two papers that are along the direction that I have in mind are:

http://onlinelibrary.wiley.com/doi/10.1002/fld.1566/abstract
http://onlinelibrary.wiley.com/doi/10.1002/fld.1723/abstract

I would love to know how to do this kind of thing on a SAMR or octree grid.

-- Boyce

Re: [petsc-users] Neumann BC with non-symmetric matrix

2016-03-01 Thread Boyce Griffith


> On Mar 1, 2016, at 12:06 PM, Mohammad Mirzadeh  wrote:
> 
> Nice discussion.
> 
> 
> On Tue, Mar 1, 2016 at 10:16 AM, Boyce Griffith  > wrote:
> 
>> On Mar 1, 2016, at 9:59 AM, Mark Adams > > wrote:
>> 
>> 
>> 
>> On Mon, Feb 29, 2016 at 5:42 PM, Boyce Griffith > > wrote:
>> 
>>> On Feb 29, 2016, at 5:36 PM, Mark Adams >> > wrote:
>>> 
>>> 
>>> GAMG is use for AMR problems like this a lot in BISICLES.
>>> 
>>> Thanks for the reference. However, a quick look at their paper suggests 
>>> they are using a finite volume discretization which should be symmetric and 
>>> avoid all the shenanigans I'm going through! 
>>> 
>>> No, they are not symmetric.  FV is even worse than vertex centered methods. 
>>>  The BCs and the C-F interfaces add non-symmetry.
>> 
>> 
>> If you use a different discretization, it is possible to make the c-f 
>> interface discretization symmetric --- but symmetry appears to come at a 
>> cost of the reduction in the formal order of accuracy in the flux along the 
>> c-f interface. I can probably dig up some code that would make it easy to 
>> compare.
>> 
>> I don't know.  Chombo/Boxlib have a stencil for C-F and do F-C with 
>> refluxing, which I do not linearize.  PETSc sums fluxes at faces directly, 
>> perhaps this IS symmetric? Toby might know.
> 
> If you are talking about solving Poisson on a composite grid, then refluxing 
> and summing up fluxes are probably the same procedure.
> 
> I am not familiar with the terminology used here. What does the refluxing 
> mean?
>  
> 
> Users of these kinds of discretizations usually want to use the conservative 
> divergence at coarse-fine interfaces, and so the main question is how to set 
> up the viscous/diffusive flux stencil at coarse-fine interfaces (or, 
> equivalently, the stencil for evaluating ghost cell values at coarse-fine 
> interfaces). It is possible to make the overall discretization symmetric if 
> you use a particular stencil for the flux computation. I think this paper 
> (http://www.ams.org/journals/mcom/1991-56-194/S0025-5718-1991-1066831-5/S0025-5718-1991-1066831-5.pdf
>  
> )
>  is one place to look. (This stuff is related to "mimetic finite difference" 
> discretizations of Poisson.) This coarse-fine interface discretization winds 
> up being symmetric (although possibly only w.r.t. a weighted inner product 
> --- I can't remember the details), but the fluxes are only first-order 
> accurate at coarse-fine interfaces.
> 
> 
> Right. I think if the discretization is conservative, i.e. discretizing div 
> of grad, and is compact, i.e. only involves neighboring cells sharing a 
> common face, then it is possible to construct symmetric discretization. An 
> example, that I have used before in other contexts, is described here: 
> http://physbam.stanford.edu/~fedkiw/papers/stanford2004-02.pdf 
> 
> 
> An interesting observation is although the fluxes are only first order 
> accurate, the final solution to the linear system exhibits super convergence, 
> i.e. second-order accurate, even in L_inf. Similar behavior is observed with 
> non-conservative, node-based finite difference discretizations. 

I don't know about that --- check out Table 1 in the paper you cite, which 
seems to indicate first-order convergence in all norms.

The symmetric discretization in the Ewing paper is only slightly more 
complicated, but will give full 2nd-order accuracy in L-1 (and maybe also L-2 
and L-infinity). One way to think about it is that you are using simple linear 
interpolation at coarse-fine interfaces (3-point interpolation in 2D, 4-point 
interpolation in 3D) using a stencil that is symmetric with respect to the 
center of the coarse grid cell.

A (discrete) Green's functions argument explains why one gets higher-order 
convergence despite localized reductions in accuracy along the coarse-fine 
interface --- it has to do with the fact that errors from individual grid 
locations do not have that large of an effect on the solution, and these c-f 
interface errors are concentrated along on a lower dimensional surface in the 
domain.

-- Boyce

Re: [petsc-users] PCFactorSetUpMatSolverPackage with SNES

2016-03-01 Thread Matthew Knepley

On Tue, Mar 1, 2016 at 12:52 PM, David Knezevic 
wrote:

> Based on KSP ex52, I use PCFactorSetUpMatSolverPackage in the process of
> setting various MUMPS ictnl options. This works fine for me when I'm
> solving linear problems.
>
> I then wanted to use PCFactorSetUpMatSolverPackage with the PC from a SNES
> object. I tried to do this with the following code (after calling
> SNESCreate, SNESSetFunction, and SNESSetJacobian):
>
> KSP snes_ksp;
> SNESGetKSP(snes, _ksp);
> PC snes_pc;
> KSPGetPC(snes_ksp, _pc);
> PCFactorSetMatSolverPackage(snes_pc, MATSOLVERMUMPS);
> PCFactorSetUpMatSolverPackage(snes_pc);
>
> However, I get a segfault on the call to PCFactorSetUpMatSolverPackage in
> this case. I was wondering what I need to do to make this work?
>
> Note that I want to set the MUMPS ictnl parameters via code rather than
> via the commandline since sometimes MUMPS fails (e.g. with error -9 due to
> a workspace size that is too small) and I need to automatically re-run the
> solve with different ictnl values when this happens.
>

That is a good reason. However, I would organize this differently. I would
still set the type from the command line.
Then later in your code, after SNES is setup correctly, you get out the PC
and reset the icntl values if you have a
failure. Its very difficult to get the setup logic correct, and its one of
the most error prone parts of PETSc. RIght now,
the most reliable way to do things is to have all the information available
up front in options.

Someday, I want to write a small DFA piece in PETSc that can encode all
this logic so that simple errors people
make can be diagnosed early with nice error messages.

  Thanks,

 Matt

> Thanks,
> David
>
>
>

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

[petsc-users] PCFactorSetUpMatSolverPackage with SNES

2016-03-01 Thread David Knezevic

Based on KSP ex52, I use PCFactorSetUpMatSolverPackage in the process of
setting various MUMPS ictnl options. This works fine for me when I'm
solving linear problems.

I then wanted to use PCFactorSetUpMatSolverPackage with the PC from a SNES
object. I tried to do this with the following code (after calling
SNESCreate, SNESSetFunction, and SNESSetJacobian):

KSP snes_ksp;
SNESGetKSP(snes, _ksp);
PC snes_pc;
KSPGetPC(snes_ksp, _pc);
PCFactorSetMatSolverPackage(snes_pc, MATSOLVERMUMPS);
PCFactorSetUpMatSolverPackage(snes_pc);

However, I get a segfault on the call to PCFactorSetUpMatSolverPackage in
this case. I was wondering what I need to do to make this work?

Note that I want to set the MUMPS ictnl parameters via code rather than via
the commandline since sometimes MUMPS fails (e.g. with error -9 due to a
workspace size that is too small) and I need to automatically re-run the
solve with different ictnl values when this happens.

Thanks,
David

Re: [petsc-users] Neumann BC with non-symmetric matrix

2016-03-01 Thread Jed Brown

Mohammad Mirzadeh  writes:

> I am not familiar with the terminology used here. What does the refluxing
> mean?

The Chombo/BoxLib family of methods evaluate fluxes between coarse grid
cells overlaying refined grids, then later visit the fine grids and
reevaluate those fluxes.  The correction needs to be propagated back to
the adjoining coarse grid cell to maintain conservation.  It's an
implementation detail that they call refluxing.

> Right. I think if the discretization is conservative, i.e. discretizing div
> of grad, and is compact, i.e. only involves neighboring cells sharing a
> common face, then it is possible to construct symmetric discretization. An
> example, that I have used before in other contexts, is described here:
> http://physbam.stanford.edu/~fedkiw/papers/stanford2004-02.pdf

It's unfortunate that this paper repeats some unfounded multigrid
slander and then basically claims to have uniform convergence using
incomplete Cholesky with CG.  In reality, incomplete Cholesky is
asymptotically no better than Jacobi.

> An interesting observation is although the fluxes are only first order
> accurate, the final solution to the linear system exhibits super
> convergence, i.e. second-order accurate, even in L_inf. 

Perhaps for aligned coefficients; definitely not for unaligned
coefficients.

signature.asc
Description: PGP signature

Re: [petsc-users] [SLEPc] record against which PETSc installation it was compiled

2016-03-01 Thread Jose E. Roman


> El 29 feb 2016, a las 17:58, Denis Davydov  escribió:
> 
> 
> 
>> On 29 Feb 2016, at 18:50, Jose E. Roman  wrote:
>> 
>> 
>>> El 29 feb 2016, a las 7:45, Denis Davydov  escribió:
>>> 
>>> Dear all,
>>> 
>>> It would be good if SLEPc would store a location of PETSc used during the 
>>> build in some 
>>> config file, e.g. `slepcconf.h`, so that this information could be 
>>> retrieved by external libraries (e.g. deal.ii) 
>>> to prevent configuring with PETSc and SLEPc while SLEPc was linking to a 
>>> different PETSc installation.
>>> See the discussion here https://github.com/dealii/dealii/issues/2167
>>> 
>>> Kind regards,
>>> Denis 
>>> 
>> 
>> I have added this:
>> https://bitbucket.org/slepc/slepc/branch/jose/configure
>> 
>> However, I am not totally convinced of this change, because PETSC_DIR is 
>> then defined both in petscconf.h and slepcconf.h, so behaviour could change 
>> depending on whether user code includes one or the other in the first place.
> 
> Thanks a lot, Jose. 
> As an alternative one could call it PETSC_DIR_IN_SLEPC or alike, so that the 
> two are different.

Ok. I changed it to SLEPC_PETSC_DIR.
Jose


> 
> Regards,
> Denis.
> 
>> 
>> I will test this further before merging into master.
>> 
>> Jose

Re: [petsc-users] Neumann BC with non-symmetric matrix

2016-03-01 Thread Mohammad Mirzadeh

Nice discussion.


On Tue, Mar 1, 2016 at 10:16 AM, Boyce Griffith 
wrote:

>
> On Mar 1, 2016, at 9:59 AM, Mark Adams  wrote:
>
>
>
> On Mon, Feb 29, 2016 at 5:42 PM, Boyce Griffith 
> wrote:
>
>>
>> On Feb 29, 2016, at 5:36 PM, Mark Adams  wrote:
>>
>>
 GAMG is use for AMR problems like this a lot in BISICLES.

>>>
>>> Thanks for the reference. However, a quick look at their paper suggests
>>> they are using a finite volume discretization which should be symmetric and
>>> avoid all the shenanigans I'm going through!
>>>
>>
>> No, they are not symmetric.  FV is even worse than vertex centered
>> methods.  The BCs and the C-F interfaces add non-symmetry.
>>
>>
>> If you use a different discretization, it is possible to make the c-f
>> interface discretization symmetric --- but symmetry appears to come at a
>> cost of the reduction in the formal order of accuracy in the flux along the
>> c-f interface. I can probably dig up some code that would make it easy to
>> compare.
>>
>
> I don't know.  Chombo/Boxlib have a stencil for C-F and do F-C with
> refluxing, which I do not linearize.  PETSc sums fluxes at faces directly,
> perhaps this IS symmetric? Toby might know.
>
>
> If you are talking about solving Poisson on a composite grid, then
> refluxing and summing up fluxes are probably the same procedure.
>

I am not familiar with the terminology used here. What does the refluxing
mean?


>
> Users of these kinds of discretizations usually want to use the
> conservative divergence at coarse-fine interfaces, and so the main question
> is how to set up the viscous/diffusive flux stencil at coarse-fine
> interfaces (or, equivalently, the stencil for evaluating ghost cell values
> at coarse-fine interfaces). It is possible to make the overall
> discretization symmetric if you use a particular stencil for the flux
> computation. I think this paper (
> http://www.ams.org/journals/mcom/1991-56-194/S0025-5718-1991-1066831-5/S0025-5718-1991-1066831-5.pdf)
> is one place to look. (This stuff is related to "mimetic finite difference"
> discretizations of Poisson.) This coarse-fine interface discretization
> winds up being symmetric (although possibly only w.r.t. a weighted inner
> product --- I can't remember the details), but the fluxes are only
> first-order accurate at coarse-fine interfaces.
>
>
Right. I think if the discretization is conservative, i.e. discretizing div
of grad, and is compact, i.e. only involves neighboring cells sharing a
common face, then it is possible to construct symmetric discretization. An
example, that I have used before in other contexts, is described here:
http://physbam.stanford.edu/~fedkiw/papers/stanford2004-02.pdf

An interesting observation is although the fluxes are only first order
accurate, the final solution to the linear system exhibits super
convergence, i.e. second-order accurate, even in L_inf. Similar behavior is
observed with non-conservative, node-based finite difference
discretizations.

-- Boyce
>

Re: [petsc-users] MatGetSize for SeqBAIJ

2016-03-01 Thread Matthew Knepley

On Tue, Mar 1, 2016 at 10:56 AM, Manav Bhatia  wrote:

> Thanks. That means I am doing something goofy in setting up my matrix.
>
> I am trying to create a matrix with block size 2, and 3000 as the number
> of block rows/columns. So, I would expect an output of 6000x6000 from the
> following code, but I get 3000x3000. Is it the sequence of my function
> calls?
>
> Thanks,
> Manav
>
>
>
> PetscErrorCode   ierr;
> Mat  mat;
> PetscInt m,n;
>
> ierr = MatCreate(mpi_comm, );
> ierr = MatSetSizes(mat, 3000, 3000, 3000, 3000);
>

MatSetSizes also takes the number of rows, not blocks.

  Matt


> ierr = MatSetType(mat, MATBAIJ);
> ierr = MatSetBlockSize(mat, 2);
> ierr = MatGetSize(mat, , );
>
> std::cout << m << "  " << n << std::endl;
>
>
>
> On Mar 1, 2016, at 10:44 AM, Hong  wrote:
>
> Manav:
>>
>>Is MatGetSize for a SeqBAIJ matrix expected to return the number of
>> block rows and columns, or the total number of rows and columns (blocks
>> rows times block size)?
>>
>
> +  m - the number of global rows
> -  n - the number of global columns
> the total number of rows and columns (blocks rows times block size).
>
> Hong
>
>
>


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

Re: [petsc-users] Investigate parallel code to improve parallelism

2016-03-01 Thread Barry Smith


  If you have access to a different cluster you can try there and see if the 
communication is any better. Likely you would get better speedup on an IBM 
BlueGene since it has a good network relative to the processing power. So best 
to run on IBM.

  Barry



> On Mar 1, 2016, at 2:03 AM, TAY wee-beng  wrote:
> 
> 
> On 29/2/2016 11:21 AM, Barry Smith wrote:
>>> On Feb 28, 2016, at 8:26 PM, TAY wee-beng  wrote:
>>> 
>>> 
>>> On 29/2/2016 9:41 AM, Barry Smith wrote:
> On Feb 28, 2016, at 7:08 PM, TAY Wee Beng  wrote:
> 
> Hi,
> 
> I've attached the files for x cells running y procs. hypre is called 
> natively I'm not sure if PETSc catches it.
   So you are directly creating hypre matrices and calling the hypre solver 
 in another piece of your code?
>>> Yes because I'm using the simple structure (struct) layout for Cartesian 
>>> grids. It's about twice as fast compared to BoomerAMG
>>Understood
>> 
>>> . I can't create PETSc matrix and use the hypre struct layout, right?
In the PETSc part of the code if you compare the 2x_y to the x_y you 
 see that doubling the problem size resulted in 2.2 as much time for the 
 KSPSolve. Most of this large increase is due to the increased time in the 
 scatter which went up to 150/54.  = 2.  but the amount of 
 data transferred only increased by 1e5/6.4e4 = 1.5625  Normally I would 
 not expect to see this behavior and would not expect such a large increase 
 in the communication time.
 
 Barry
 
 
 
>>> So ideally it should be 2 instead of 2.2, is that so?
>>   Ideally
>> 
>>> May I know where are you looking at? Because I can't find the nos.
>>   The column labeled Avg len tells the average length of messages which 
>> increases from 6.4e4 to 1e5 while the time max increase by 2.77 (I took the 
>> sum of the VecScatterBegin and VecScatter End rows.
>> 
>>> So where do you think the error comes from?
>>   It is not really an error it is just that it is taking more time then one 
>> would hope it would take.
>>> Or how can I troubleshoot further?
>> 
>>If you run the same problem several times how much different are the 
>> numerical timings for each run?
> Hi,
> 
> I have re-done x_y and 2x_y again. I have attached the files with _2 for the 
> 2nd run. They're exactly the same.
> 
> Should I try running on another cluster?
> 
> I also tried running the same problem with more cells and more time steps (to 
> reduce start up effects) on another cluster. But I forgot to run it with 
> -log_summary. Anyway, the results show:
> 
> 1. Using 1.5 million cells with 48 procs and 3M with 96p took 65min and 
> 69min. Using the weak scaling formula I attached earlier, it gives about 88% 
> efficiency
> 
> 2. Using 3 million cells with 48 procs and 6M with 96p took 114min and 
> 121min. Using the weak scaling formula I attached earlier, it gives about 88% 
> efficiency
> 
> 3. Using 3.75 million cells with 48 procs and 7.5M with 96p took 134min and 
> 143min. Using the weak scaling formula I attached earlier, it gives about 87% 
> efficiency
> 
> 4. Using 4.5 million cells with 48 procs and 9M with 96p took 160min and 
> 176min (extrapolated). Using the weak scaling formula I attached earlier, it 
> gives about 80% efficiency
> 
> So it seems that I should run with 3.75 million cells with 48 procs and scale 
> along this ratio. Beyond that, my efficiency decreases. Is that so? Maybe I 
> should also run with -log_summary to get better estimate...
> 
> Thanks.
>> 
>> 
>>> Thanks
> Thanks
> 
> On 29/2/2016 1:11 AM, Barry Smith wrote:
>>   As I said before, send the -log_summary output for the two processor 
>> sizes and we'll look at where it is spending its time and how it could 
>> possibly be improved.
>> 
>>   Barry
>> 
>>> On Feb 28, 2016, at 10:29 AM, TAY wee-beng  wrote:
>>> 
>>> 
>>> On 27/2/2016 12:53 AM, Barry Smith wrote:
> On Feb 26, 2016, at 10:27 AM, TAY wee-beng  wrote:
> 
> 
> On 26/2/2016 11:32 PM, Barry Smith wrote:
>>> On Feb 26, 2016, at 9:28 AM, TAY wee-beng  wrote:
>>> 
>>> Hi,
>>> 
>>> I have got a 3D code. When I ran with 48 procs and 11 million 
>>> cells, it runs for 83 min. When I ran with 96 procs and 22 million 
>>> cells, it ran for 99 min.
>>This is actually pretty good!
> But if I'm not wrong, if I increase the no. of cells, the parallelism 
> will keep on decreasing. I hope it scales up to maybe 300 - 400 procs.
>>> Hi,
>>> 
>>> I think I may have mentioned this before, that is, I need to submit a 
>>> proposal to request for computing nodes. In the proposal, I'm supposed 
>>> to run some simulations to estimate the time it takes to run my code.

Re: [petsc-users] MatGetSize for SeqBAIJ

2016-03-01 Thread Manav Bhatia

Thanks. That means I am doing something goofy in setting up my matrix. 

I am trying to create a matrix with block size 2, and 3000 as the number of 
block rows/columns. So, I would expect an output of 6000x6000 from the 
following code, but I get 3000x3000. Is it the sequence of my function calls? 

Thanks,
Manav

PetscErrorCode   ierr;
Mat  mat;
PetscInt m,n;

ierr = MatCreate(mpi_comm, );
ierr = MatSetSizes(mat, 3000, 3000, 3000, 3000); 
ierr = MatSetType(mat, MATBAIJ); 
ierr = MatSetBlockSize(mat, 2);  
ierr = MatGetSize(mat, , );

std::cout << m << "  " << n << std::endl;

> On Mar 1, 2016, at 10:44 AM, Hong  wrote:
> 
> Manav:
>Is MatGetSize for a SeqBAIJ matrix expected to return the number of block 
> rows and columns, or the total number of rows and columns (blocks rows times 
> block size)?
>  
> +  m - the number of global rows
> -  n - the number of global columns
> the total number of rows and columns (blocks rows times block size).
> 
> Hong
>

Re: [petsc-users] MatGetSize for SeqBAIJ

2016-03-01 Thread Hong

Manav:
>
>Is MatGetSize for a SeqBAIJ matrix expected to return the number of
> block rows and columns, or the total number of rows and columns (blocks
> rows times block size)?
>

+  m - the number of global rows
-  n - the number of global columns
the total number of rows and columns (blocks rows times block size).

Hong

[petsc-users] MatGetSize for SeqBAIJ

2016-03-01 Thread Manav Bhatia

Hi

   Is MatGetSize for a SeqBAIJ matrix expected to return the number of block 
rows and columns, or the total number of rows and columns (blocks rows times 
block size)?

Thanks,
Manav

Re: [petsc-users] how to get full trajectory from TS into petsc binary

2016-03-01 Thread Barry Smith


> On Mar 1, 2016, at 12:39 AM, Ed Bueler  wrote:
> 
> Barry --
> 
> Will try it.
> 
> > ... since, presumably, other more powerful IO tools exist that would be 
> > used for "real" problems?
> 
> I know there are tools for snapshotting from PETSc, e.g. VecView to .vtk.  In 
> fact petsc binary seems fairly convenient for that.  On the other hand, I am 
> not sure I've ever done anything "real".  ;-)
> 
> Anyone out there:  Are there a good *convenient* tools for saving 
> space/time-series (= movies) from PETSc TS?  I want to add frames and movies 
> from PETSc into slides, etc.  I can think of NetCDF but it seems 
> not-very-convenient, and I am worried not well-supported from PETSc.  Is 
> setting up TS with events (=TSSetEventMonitor()) and writing separate 
> snapshot files the preferred scalable usage, despite the extra effort 
> compared to "-ts_monitor_solution binary:foo.dat"?

   Ed,

As I said in my previous email since we don't have a way of indicating 
plain double precision numbers in our binary files it is not possible to put 
both the vectors and the time steps in the same file without augmenting our 
file format.

  Barry



> 
> Ed
> 
> 
> On Mon, Feb 29, 2016 at 8:53 PM, Barry Smith  wrote:
> 
>   Ed,
> 
>I have added a branch barry/feature-ts-monitor-binary  that supports 
> -ts_monitor binary:timesteps that will store in simple binary format each of 
> the time steps associated with each solution. This in conjugation with 
> -ts_monitor_solution binary:solutions will give you two files you can read 
> in. But note that timesteps is a simple binary file of double precision 
> numbers you should read in directly in python, you cannot use 
> PetscBinaryIO.py which is what you will use to read in the solutions file.
> 
>   Barry
> 
> Currently PETSc has a binary file format where we can save Vec, Mat, IS, each 
> is marked with a type id for PetscBinaryIO.py to detect, we do not have type 
> ids for simple double precision numbers or arrays of numbers. This is why I 
> have no way of saving the time steps in a way that PetscBinaryIO.py could 
> read them in currently. I don't know how far we want to go in "spiffing up" 
> the PETSc binary format to do more elaborate things since, presumably, other 
> more power IO tools exist that would be used for "real" problems?
> 
> 
> > On Feb 29, 2016, at 3:24 PM, Ed Bueler  wrote:
> >
> > Dear PETSc --
> >
> > I have a short C ode code that uses TS to solve  y' = g(t,y)  where y(t) is 
> > a 2-dim'l vector.  My code defaults to -ts_type rk so it does adaptive 
> > time-stepping; thus using -ts_monitor shows times at stdout:
> >
> > $ ./ode -ts_monitor
> > solving from t0 = 0.000 with initial time step dt = 0.1 ...
> > 0 TS dt 0.1 time 0.
> > 1 TS dt 0.170141 time 0.1
> > 2 TS dt 0.169917 time 0.270141
> > 3 TS dt 0.171145 time 0.440058
> > 4 TS dt 0.173931 time 0.611203
> > 5 TS dt 0.178719 time 0.785134
> > 6 TS dt 0.0361473 time 0.963853
> > 7 TS dt 0.188252 time 1.
> > error at tf = 1.000 :  |y-y_exact|_inf = 0.000144484
> >
> > I want to output the trajectory in PETSc binary and plot it in python using 
> > bin/PetscBinaryIO.py.  Clearly I need the times shown above to do that.
> >
> > Note "-ts_monitor_solution binary:XX" gives me a binary file with only y 
> > values in it, but not the corresponding times.
> >
> > My question is, how to get those times in either the same binary file 
> > (preferred) or separate binary files?  I have tried
> >
> > $ ./ode -ts_monitor binary:foo.dat# invalid
> > $ ./ode -ts_monitor_solution binary:bar.dat   # no t in file
> > $ ./ode -ts_monitor_solution binary:baz.dat -ts_save_trajectory   # no t in 
> > file
> >
> > without success.  (I am not sure what the boolean option 
> > -ts_save_trajectory does, by the way.)
> >
> > Thanks!
> >
> > Ed
> >
> > PS Sorry if this is a "RTFM" question, but so far I can't find the 
> > documentation.
> >
> >
> > --
> > Ed Bueler
> > Dept of Math and Stat and Geophysical Institute
> > University of Alaska Fairbanks
> > Fairbanks, AK 99775-6660
> > 301C Chapman and 410D Elvey
> > 907 474-7693 and 907 474-7199  (fax 907 474-5394)
> 
> 
> 
> 
> -- 
> Ed Bueler
> Dept of Math and Stat and Geophysical Institute
> University of Alaska Fairbanks
> Fairbanks, AK 99775-6660
> 301C Chapman and 410D Elvey
> 907 474-7693 and 907 474-7199  (fax 907 474-5394)

Re: [petsc-users] Read and use boundaries from exodus in finite element context

2016-03-01 Thread Matthew Knepley

On Tue, Mar 1, 2016 at 9:58 AM, Max Rietmann  wrote:

> Dear all,
>
> We are building out a new finite-element code for wave propagation and I
> am currently implementing the boundary conditions. My first pass will
> simply have dirichlet boundaries, but eventually we will have more
> sophisticated options available.
>

Is this going to be a purely explicit code?


> I am creating an exodus mesh in Trelis/Cubit, in which I can create one
> (or more) "side sets", which are properly read by the exodus reader. From
> reading the petsc source (plexexodusii.c), it seems that these basically
> create a list of faces, which belong to the side set.
>

Yes.


> A call to DMPlexGetLabelValue(dm,"Face Set",face_id,), allows me to
> see if "face_id" is on the boundary by checking value (>=1 for boundary, or
> -1 for not in set). Additional side sets get ids = {2,3,etc}. This allows
> us to have multiple types of boundary (absorbing, reflecting, etc).
>

Yes. Note that you can also extract them all at once into a sorted array.


> However, I am unsure how to determine if an element has a particular
> face_id in order to determine if one face of the element is on a boundary
> (and which value it has {-1,1,2,3, etc}).
>

Here you can see me doing the thing you are asking for:


https://bitbucket.org/petsc/petsc/src/2d7a10a4145949cccd1c2cb7dc0d518dc12666a9/src/dm/impls/plex/plexsubmesh.c?at=master=file-view-default#plexsubmesh.c-106

The code is slightly overkill for what you want since it also works for
vertices. You could call

  DMPlexGetSupport()

on the face rather than

  DMPlexGetTransitiveClosure()

and the code would be simpler.


> The routine is listed here:
>
> PetscErrorCode DMPlexGetLabelValue(DM dm, const char name[], PetscInt
> point, PetscInt *value)
>
> How do I determine the "point" of a face, if I'm iterating through my
> elements.
>
> thanks!
>
> Max Rietmann
>
> PS I've also seen the DMPlexAddBoundary method, but I wasn't sure how to
> use it in our setting.
>

This is used to hook into the DMPlexSetBoundaryValues() interface so that
BC values are properly loaded into
local vectors before assembly operations.

  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

[petsc-users] Read and use boundaries from exodus in finite element context

2016-03-01 Thread Max Rietmann


Dear all,

We are building out a new finite-element code for wave propagation and I 
am currently implementing the boundary conditions. My first pass will 
simply have dirichlet boundaries, but eventually we will have more 
sophisticated options available.


I am creating an exodus mesh in Trelis/Cubit, in which I can create one 
(or more) "side sets", which are properly read by the exodus reader. 
From reading the petsc source (plexexodusii.c), it seems that these 
basically create a list of faces, which belong to the side set.


A call to DMPlexGetLabelValue(dm,"Face Set",face_id,), allows me 
to see if "face_id" is on the boundary by checking value (>=1 for 
boundary, or -1 for not in set). Additional side sets get ids = 
{2,3,etc}. This allows us to have multiple types of boundary (absorbing, 
reflecting, etc).


However, I am unsure how to determine if an element has a particular 
face_id in order to determine if one face of the element is on a 
boundary (and which value it has {-1,1,2,3, etc}).


The routine is listed here:

PetscErrorCode DMPlexGetLabelValue(DM dm, const char name[], 
PetscInt point, PetscInt *value)


How do I determine the "point" of a face, if I'm iterating through my 
elements.


thanks!

Max Rietmann

PS I've also seen the DMPlexAddBoundary method, but I wasn't sure how to 
use it in our setting.

Re: [petsc-users] Neumann BC with non-symmetric matrix

2016-03-01 Thread Jed Brown

Boyce Griffith  writes:
> Jed, can you also do this for Stokes?  It seems like something like
> RT0 is the right place to start.

See, for example, Arnold, Falk, and Winther's 2007 paper on mixed FEM
for elasticity with weakly imposed symmetry.  It's the usual H(div)
methodology and should apply equally well to Stokes.  I'm not aware of
any analysis or results of choosing quadrature to eliminate flux terms
in these discretizations.

signature.asc
Description: PGP signature

Re: [petsc-users] Neumann BC with non-symmetric matrix

2016-03-01 Thread Boyce Griffith


> On Mar 1, 2016, at 10:08 AM, Jed Brown  wrote:
> 
> Mark Adams  writes:
> 
>> I don't know.  Chombo/Boxlib have a stencil for C-F and do F-C with
>> refluxing, which I do not linearize.  PETSc sums fluxes at faces directly,
>> perhaps this IS symmetric? 
> 
> Depends on scaling.  Refluxing versus evaluating fluxes once on the fine
> faces should give the same result.  Note that you can make a symmetric
> FV discretization (for elliptic problems) by writing it as a mixed FE
> method and choosing a quadrature so that fluxes can be eliminated.

Jed, can you also do this for Stokes?  It seems like something like RT0 is the 
right place to start.

Thanks,

-- Boyce

Re: [petsc-users] Example for FAS

2016-03-01 Thread Cyrill Vonplanta

Thanks a lot. Checking out with example 5 works fine, problem solved (for now).

Cyrill

Re: [petsc-users] Neumann BC with non-symmetric matrix

2016-03-01 Thread Boyce Griffith

> On Mar 1, 2016, at 9:59 AM, Mark Adams  wrote:
> 
> 
> 
> On Mon, Feb 29, 2016 at 5:42 PM, Boyce Griffith  > wrote:
> 
>> On Feb 29, 2016, at 5:36 PM, Mark Adams > > wrote:
>> 
>> 
>> GAMG is use for AMR problems like this a lot in BISICLES.
>> 
>> Thanks for the reference. However, a quick look at their paper suggests they 
>> are using a finite volume discretization which should be symmetric and avoid 
>> all the shenanigans I'm going through!
>> 
>> No, they are not symmetric.  FV is even worse than vertex centered methods.  
>> The BCs and the C-F interfaces add non-symmetry.
> 
> 
> If you use a different discretization, it is possible to make the c-f 
> interface discretization symmetric --- but symmetry appears to come at a cost 
> of the reduction in the formal order of accuracy in the flux along the c-f 
> interface. I can probably dig up some code that would make it easy to compare.
> 
> I don't know.  Chombo/Boxlib have a stencil for C-F and do F-C with 
> refluxing, which I do not linearize.  PETSc sums fluxes at faces directly, 
> perhaps this IS symmetric? Toby might know.

If you are talking about solving Poisson on a composite grid, then refluxing 
and summing up fluxes are probably the same procedure.

Users of these kinds of discretizations usually want to use the conservative 
divergence at coarse-fine interfaces, and so the main question is how to set up 
the viscous/diffusive flux stencil at coarse-fine interfaces (or, equivalently, 
the stencil for evaluating ghost cell values at coarse-fine interfaces). It is 
possible to make the overall discretization symmetric if you use a particular 
stencil for the flux computation. I think this paper 
(http://www.ams.org/journals/mcom/1991-56-194/S0025-5718-1991-1066831-5/S0025-5718-1991-1066831-5.pdf

)
 is one place to look. (This stuff is related to "mimetic finite difference" 
discretizations of Poisson.) This coarse-fine interface discretization winds up 
being symmetric (although possibly only w.r.t. a weighted inner product --- I 
can't remember the details), but the fluxes are only first-order accurate at 
coarse-fine interfaces.

-- Boyce

Re: [petsc-users] Neumann BC with non-symmetric matrix

2016-03-01 Thread Jed Brown

Mark Adams  writes:

> I don't know.  Chombo/Boxlib have a stencil for C-F and do F-C with
> refluxing, which I do not linearize.  PETSc sums fluxes at faces directly,
> perhaps this IS symmetric? 

Depends on scaling.  Refluxing versus evaluating fluxes once on the fine
faces should give the same result.  Note that you can make a symmetric
FV discretization (for elliptic problems) by writing it as a mixed FE
method and choosing a quadrature so that fluxes can be eliminated.

signature.asc
Description: PGP signature

Re: [petsc-users] Neumann BC with non-symmetric matrix

2016-03-01 Thread Mark Adams

On Mon, Feb 29, 2016 at 5:42 PM, Boyce Griffith 
wrote:

>
> On Feb 29, 2016, at 5:36 PM, Mark Adams  wrote:
>
>
>>> GAMG is use for AMR problems like this a lot in BISICLES.
>>>
>>
>> Thanks for the reference. However, a quick look at their paper suggests
>> they are using a finite volume discretization which should be symmetric and
>> avoid all the shenanigans I'm going through!
>>
>
> No, they are not symmetric.  FV is even worse than vertex centered
> methods.  The BCs and the C-F interfaces add non-symmetry.
>
>
> If you use a different discretization, it is possible to make the c-f
> interface discretization symmetric --- but symmetry appears to come at a
> cost of the reduction in the formal order of accuracy in the flux along the
> c-f interface. I can probably dig up some code that would make it easy to
> compare.
>

I don't know.  Chombo/Boxlib have a stencil for C-F and do F-C with
refluxing, which I do not linearize.  PETSc sums fluxes at faces directly,
perhaps this IS symmetric? Toby might know.


>
> -- Boyce
>

Re: [petsc-users] Neumann BC with non-symmetric matrix

2016-03-01 Thread Jed Brown

Barry Smith  writes:
>>4-2) Another observation is that the true residual stagnates way
>>earlier which I assume is an indication that the RHS is in fact
>>not in the range of A.
>
>You can hope this is the case.
>
>Of course one cannot really know if the residual is stagnating due
>to an inconsistent right hand side or for some other reason. But if
>it stagnates at 10^-4 it is probably due to inconsistent right hand
>side if it stagnates at 10^-12 the right hand side is probably
>consistent. If it stagnates at 10^-8 then it could be due to
>inconsistent rhs and or some other reason.

If the true residual stagnates while the preconditioned residual
continues to converge, it may be that the preconditioner is singular or
nearly so.  This often happens if you use a black-box method for a
saddle point problem, for example.


signature.asc
Description: PGP signature

Re: [petsc-users] Example for FAS

2016-03-01 Thread Jed Brown

Cyrill Vonplanta  writes:

> Dear PETSc User-list,
>
> I am trying to use the FAS from PETSc in version 3.6.3. For other solvers 
> there would be usually an example or two on the documentation, but on 
> http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/SNES/SNESFAS.html 
> I couldn’t find one at the bottom.

It's pure run-time options.  See src/snes/examples/tutorials/makefile
and search for '-snes_type fas'.  ex5 and ex19 are good examples to
check out.

signature.asc
Description: PGP signature

Re: [petsc-users] how to get full trajectory from TS into petsc binary

2016-03-01 Thread Jed Brown

Ed Bueler  writes:

> Barry --
>
> I am reading the resulting file successfully using
>
> import struct
> import sys
> f = open('timesteps','r')
> while True:
> try:
> bytes = f.read(8)
> except:
> print "f.read() failed"
> sys.exit(1)

It seems odd to catch the exception and convert to sys.exit(1).  This
just loses the stack trace (which might have more detailed information)
and makes the debugger catch the wrong location.

> if len(bytes) > 0:
> print struct.unpack('>d',bytes)[0]
> else:
> break
> f.close()
>
> However, was there a more elegant intended method?  

I would use numpy.fromfile, which is one line.

> I am disturbed by the apparent need to specify big-endian-ness (=
> '>d') in the struct.unpack() call.

You need it because PETSc binary files are big-endian.  There is no
place to encode the byte order in these raw binary files, so some
convention is required for the file to portable.

HDF5 includes this information, though it is an annoying dependency.
Numpy files (*.npy) are a simpler alternative that PETSc could decide to
support.

https://docs.scipy.org/doc/numpy/neps/npy-format.html

signature.asc
Description: PGP signature

[petsc-users] Example for FAS

2016-03-01 Thread Cyrill Vonplanta

Dear PETSc User-list,

I am trying to use the FAS from PETSc in version 3.6.3. For other solvers there 
would be usually an example or two on the documentation, but on 
http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/SNES/SNESFAS.html I 
couldn’t find one at the bottom.

Does anybody have an example for FAS using different preconditioners etc.?

Best regards
Cyrill

Re: [petsc-users] MOOSE_SNES_ISSUES

2016-03-01 Thread Patrick Sanan

On Tue, Mar 01, 2016 at 09:58:18AM +0100, alena kopanicakova wrote:
> Hello,
> 
> I am developing my own nonlinear solver, which aims to solve simulations
> from MOOSE. Communication with moose is done via SNES interface.
> 
> I am obtaining Jacobian and residual in following way:
> 
>SNESComputeFunction(snes, X, R);
> 
>   SNESSetJacobian(snes, jac, jac, SNESComputeJacobianDefault,
As here, 
http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/SNES/SNESComputeJacobianDefault.html
 , 
it seems like you are computing the Jacobian in a brute-force way with finite 
differences. 
Do you have an expression to compute the Jacobian?
> NULL);
>   SNESComputeJacobian(snes, X, jac,  jac);
> 
> Unfortunately, by this setting it takes incredible amount of time to obtain
> Jacobian.
> Taking closer look at perf log, it seems that difference between mine and
> MOOSE solver is, that my executioner calls compute_residual() function
> ridiculously many times.
> I have no idea, what could be causing such a behavior.
> 
> Do you have any suggestions, how to set up interface properly? or which
> change is needed for obtaining more-less same performance as MOOSE
> executioner?
> 
> 
> many thanks,
> alena

> 
> Framework Information:
> MOOSE version:   git commit e79292e on 2016-01-05
> PETSc Version:   3.6.1
> Current Time:Tue Mar  1 00:51:45 2016
> Executable Timestamp:Tue Mar  1 00:36:30 2016
> 
> Parallelism:
>   Num Processors:  1
>   Num Threads: 1
> 
> Mesh: 
>   Distribution:serial
>   Mesh Dimension:  3
>   Spatial Dimension:   3
>   Nodes:   
> Total: 1331
> Local: 1331
>   Elems:   
> Total: 1000
> Local: 1000
>   Num Subdomains:  1
>   Num Partitions:  1
> 
> Nonlinear System:
>   Num DOFs:3993
>   Num Local DOFs:  3993
>   Variables:   { "disp_x" "disp_y" "disp_z" } 
>   Finite Element Types:"LAGRANGE" 
>   Approximation Orders:"FIRST" 
> 
> Execution Information:
>   Executioner: Steady
>   Solver Mode: NEWTON
> 
> 
> 
>  0 Nonlinear |R| = [32m2.85e-02[39m
>   0 Linear |R| = [32m2.85e-02[39m
>   1 Linear |R| = [32m1.653670e-02[39m
>   2 Linear |R| = [32m1.447168e-02[39m
>   3 Linear |R| = [33m1.392965e-02[39m
>   4 Linear |R| = [32m1.258440e-02[39m
>   5 Linear |R| = [32m1.007181e-02[39m
>   6 Linear |R| = [32m8.264315e-03[39m
>   7 Linear |R| = [32m6.541897e-03[39m
>   8 Linear |R| = [32m4.371900e-03[39m
>   9 Linear |R| = [32m2.100406e-03[39m
>  10 Linear |R| = [32m1.227539e-03[39m
>  11 Linear |R| = [32m1.026286e-03[39m
>  12 Linear |R| = [32m9.180101e-04[39m
>  13 Linear |R| = [32m8.087598e-04[39m
>  14 Linear |R| = [32m6.435247e-04[39m
>  15 Linear |R| = [32m5.358688e-04[39m
>  16 Linear |R| = [32m4.551657e-04[39m
>  17 Linear |R| = [32m4.090276e-04[39m
>  18 Linear |R| = [32m3.359290e-04[39m
>  19 Linear |R| = [32m2.417643e-04[39m
>  20 Linear |R| = [32m1.710452e-04[39m
>  21 Linear |R| = [32m1.261996e-04[39m
>  22 Linear |R| = [32m9.384052e-05[39m
>  23 Linear |R| = [32m6.070637e-05[39m
>  24 Linear |R| = [32m4.283233e-05[39m
>  25 Linear |R| = [32m3.383792e-05[39m
>  26 Linear |R| = [32m2.342289e-05[39m
>  27 Linear |R| = [32m1.700855e-05[39m
>  28 Linear |R| = [32m9.814278e-06[39m
>  29 Linear |R| = [32m4.398519e-06[39m
>  30 Linear |R| = [32m2.161205e-06[39m
>  31 Linear |R| = [32m1.289206e-06[39m
>  32 Linear |R| = [32m6.548007e-07[39m
>  33 Linear |R| = [32m3.677894e-07[39m
>  34 Linear |R| = [32m2.640006e-07[39m
>  1 Nonlinear |R| = [32m2.400310e-02[39m
>   0 Linear |R| = [32m2.400310e-02[39m
>   1 Linear |R| = [32m9.102075e-03[39m
>   2 Linear |R| = [32m5.017381e-03[39m
>   3 Linear |R| = [32m3.840732e-03[39m
>   4 Linear |R| = [32m2.990685e-03[39m
>   5 Linear |R| = [32m1.990203e-03[39m
>   6 Linear |R| = [32m1.085764e-03[39m
>   7 Linear |R| = [32m4.657779e-04[39m
>   8 Linear |R| = [32m3.049692e-04[39m
>   9 Linear |R| = [32m1.625839e-04[39m
>  10 Linear |R| = [32m1.124700e-04[39m
>  11 Linear |R| = [32m7.764153e-05[39m
>  12 Linear |R| = [32m5.698577e-05[39m
>  13 Linear |R| = [32m4.581843e-05[39m
>  14 Linear |R| = [32m4.262610e-05[39m
>  15 Linear |R| = [32m3.792804e-05[39m
>  16 Linear |R| = [32m3.404168e-05[39m
>  17 Linear |R| = [32m2.536004e-05[39m
>  18 Linear |R| = [32m1.577559e-05[39m
>  19 Linear |R| = [32m9.099392e-06[39m
>  20 Linear |R| = [32m6.140685e-06[39m
>  21 Linear |R| = [32m5.083606e-06[39m
>  22 Linear |R| =

37 matches

Mail list logo