[petsc-users] Tracking a NaN in FVM gradient computation

[Ellen M. Price]

Hi there PETSc,

I am working my way through ex11.c and have encountered a problem. On
the first pass through mesh adaption, the gradient computation in PETSc
triggers a NaN, even though none of the input data are NaN.

To reproduce:

PETSc v3.15 with latest MPICH, no external libraries, debugging on
Compile ex11.c using GCC 9.3 on Ubuntu 20.04
Run as: ./ex11 -ufv_use_amr -fp_trap

Running under GDB shows that the offending line is:
src/dm/impls/plex/plexfvm.c:111
This originates from the call to DMPlexReconstructGradientsFVM.

I only noticed this because I was trying to resolve a discontinuity in
my own initial condition for a Sod shock tube problem, but then I found
that it occurs in the example without any modification, too.

Is this somehow intended? If not, what steps can one take to make sure
the gradients are actual numeric values? It doesn't make sense to me
that they would be undefined on the first step when nothing else is NaN.
Or is it that AMR in ex11 doesn't quite work yet?

Thanks,
Ellen Price

[petsc-users] Guidance on using Tao for unconstrained minimization

[Ellen M. Price]

Hello PETSc users!

I am using Tao for an unconstrained minimization problem. I have found
that CG works better than the other types for this application. After
about 85 iterations, I get an error about line search failure. I'm not
clear on what this means, or how I could mitigate the problem, and
neither the manual nor FAQ give any guidance. Can anyone suggest things
I could try to help the method converge? I have function and gradient
info, but no Hessian.

Thanks,
Ellen Price

Re: [petsc-users] Updating TS solution outside PETSc

[Ellen M. Price]

I think I'm still unclear on exactly how this should work. Suppose, in
my RHS function for TS, I'm processing the grid to compute its time
derivative and get to an edge. What do I do?

If I set the derivative there to zero, the value will never change, but
it *should* change so that the spatial derivative there is zero.

If I set it to the value it would get if it wasn't an edge, then the
derivative isn't preserved anymore.

This is where I get stuck.

Ellen


On 12/5/19 10:16 AM, Smith, Barry F. wrote:
> 
>Are you using cell-centered or vertex centered discretization ( makes a 
> slight difference)?
> 
>Our model is to use DM_BOUNDARY_MIRROR  DMBoundaryType. This means that 
> u_first_real_grid_point - u_its_ghost_point = 0 (since DMGlobalToLocal will 
> automatically put into the physical ghost location the appropriate mirror 
> values) thus u_n is zero along the edge; zero Neumann  conditions, for 
> non-zero Neuman you need to put something in the "local rhs" to represent 
> that, I'm not sure it is clear, but think about it in one dimension for the 
> non-zero Neumann case.
> 
>Bad news  not yet implemented for 3d.
> 
>If you are using 3d we should fix this for you (or you fix it and make a 
> MR because we should have this support). 
> 
>If you have a 2d version of your code I would test with 3d your model etc 
> and then let us know and request how we can get the 3d mirror written.
> 
>Others may have alternative advice for Neumann with DMDA,
> 
>Barry
> 
> 
>> On Dec 5, 2019, at 10:00 AM, Ellen M. Price  
>> wrote:
>>
>> Hi PETSc users,
>>
>> I am working with a code that solves a set of PDEs on a rectangular
>> domain with Neumann boundary conditions. My understanding of
>> implementing the boundary condition is that I should set the boundary
>> value to be the value that makes the finite difference derivative go to
>> zero (or some other prescribed value) on that boundary.
>>
>> I was attempting to update the solution from TS using a pre- or
>> post-step/stage function and TSGetSolution, but this does not appear to
>> work as expected. What would be the correct way to prescribe the
>> boundary condition, given that I'm using TSRK for timestepping and a
>> DMDA for discretization?
>>
>> Looking forward to any help!
>>
>> Ellen Price
> 

[petsc-users] Updating TS solution outside PETSc

[Ellen M. Price]

Hi PETSc users,

I am working with a code that solves a set of PDEs on a rectangular
domain with Neumann boundary conditions. My understanding of
implementing the boundary condition is that I should set the boundary
value to be the value that makes the finite difference derivative go to
zero (or some other prescribed value) on that boundary.

I was attempting to update the solution from TS using a pre- or
post-step/stage function and TSGetSolution, but this does not appear to
work as expected. What would be the correct way to prescribe the
boundary condition, given that I'm using TSRK for timestepping and a
DMDA for discretization?

Looking forward to any help!

Ellen Price

[petsc-users] DMCompositeSetCoupling -- simple example?

[Ellen M. Price via petsc-users]

Hello again!

For my multiphysics problem, I think a DMComposite might make the most
sense, since some variables are only defined in 1d (so any 2d
information is just wasting memory and CPU), while others have a full 2d
distribution. To tell PETSc how the two DMs should be coupled, it looks
like I need DMCompositeSetCoupling, but there are no examples on how to
use it, or even what the parameters of the provided function mean. Does
anyone have a simple example to share? Doesn't have to be full code, a
snippet would be just fine.

Thanks,
Ellen Price

[petsc-users] Correct way to access a sequential version of a DMDA Vec?

[Ellen M. Price via petsc-users]

Hi PETSc users!

I'm trying to access a DMDA Vec's values after doing a scatter to all
processors. I've run into some trouble, however, because the call to
DMDAVecGetArray doesn't seem to care that I'm using a sequential vector.
The calling sequence that "works" (runs without error *until* I try to
use the values) is:

// da is the 2-dimensional DMDA
// x is a global, input vector
// xs is the sequential vector
// values is an array of structs suitable for this DMDA

DMGetLocalVector(da, &xloc);
DMGlobalToLocalBegin(da, x, INSERT_VALUES, xloc);
DMGlobalToLocalEnd(da, x, INSERT_VALUES, xloc);
VecScatterCreateToAll(xloc, &scatter, &xs);
VecScatterBegin(scatter, xloc, xs, INSERT_VALUES, SCATTER_FORWARD);
VecScatterEnd(scatter, xloc, xs, INSERT_VALUES, SCATTER_FORWARD);
DMDAVecGetArrayRead(da, xs, &values);

>From a skim of the PETSc source code, it looks like the start parameters
that get handed off to VecGetArray2d only depend on the DMDA and do not
take into account that the vector might already be sequential. (I might
be wrong about this, but that's what it looks like at first glance.) Is
there a way around this that is more elegant than just modifying the
output pointer?

Feel free to point me to an applicable message board post if there is
one, but Google hasn't brought up anything useful yet.

Thanks,
Ellen Price

Re: [petsc-users] Changing nonzero structure and Jacobian coloring

[Ellen M. Price via petsc-users]

 the problem using the background search cell grid and do not 
>> load balance the particles which is a commonly done with incompressible SPH 
>> implementations). As a result, the local size of the Vec object which stores 
>> the solution will change between time steps. Vec cannot be re-sized, hence 
>> you will not only need to change the nonzero structure of the Jacobian but 
>> you will also need to destroy/create all Vec's objects and all objects 
>> associated with the nonlinear solve. Given this, I'm not even sure you can 
>> use TS for your use case (hopefully a TS expert will comment on this). 
>>
>> My experience has been that creating new objects (Vec, Mat, KSP, SNES) in 
>> PETSc is fast (compared to a linear solve). You might have to give up on 
>> using TS, and instead roll your own time integrator. By doing this you will 
>> have control of only a SNES object (plus Mat's for J Vec's for the residual 
>> and solution) which you can create and destroy within each time step. To use 
>> FD coloring you would need to provide this function 
>> SNESComputeJacobianDefaultColor to SNESComputeJacobian(), along with a 
>> preallocated matrix (which you define using MatPreallocator).
>>
>> If rolling your own TS would work, then TSStep() will definitely work, and 
>> saves a lot of coding for multistep methods.
>>
>>   Thanks,
>>
>> Matt
>>  
>> Thanks
>> Dave
>> Dave
>>
>>
>>
>> On Wed 16. Oct 2019 at 13:25, Matthew Knepley via petsc-users 
>>  wrote:
>> On Tue, Oct 15, 2019 at 4:56 PM Smith, Barry F. via petsc-users 
>>  wrote:
>>
>>   Because of the difficulty of maintaining a nonzero matrix for such 
>> problems aren't these problems often done "matrix-free"?
>>
>>   So you provide a routine that computes the action of the operator but 
>> doesn't form the operator explicitly (and you hope that you don't need a 
>> preconditioner). There are simple ways (but less optimal) to do this as well 
>> as more sophisticated (such as multipole methods). 
>>
>>   If the convergence of the linear solver is too slow (due to lack of 
>> preconditioner) you might consider continuing with matrix free but at each 
>> new Newton solve (or several Newton solves) construct a very sparse matrix 
>> that captures just the very local coupling in the problem. Once particles 
>> have moved around a bit you would throw away the old matrix and construct a 
>> new one. For example the matrix might just captures interactions between 
>> particles that are less than some radius away from themselves. You could use 
>> a direct solver or iterative solver to solve this very sparse system.
>>
>> I tried to do this with Dan Negrut many years ago and had the same problems. 
>> That is part of why incompressible SPH never works,
>> since you need global modes.
>>
>>   Thanks,
>>
>>  Matt
>>  
>>   Barry
>>
>> This is why KSPSetOperators and SNESSetJacobian and TSSetRHS/IJacobain take 
>> two Jacobian matrices, the first holds the matrix free Jacobian and the 
>> second holds the approximation used inside the preconditioner.
>>
>>
>>
>>> On Oct 15, 2019, at 3:29 PM, Ellen M. Price  
>>> wrote:
>>>
>>> Thanks for the reply, Barry! Unfortunately, this is a particle code
>>> (SPH, specifically), so the particle neighbors, which influence the
>>> properties, change over time; the degrees of freedom is a constant, as
>>> is the particle number. Any thoughts, given the new info? Or should I
>>> stick with explicit integration? I've seen explicit used most commonly,
>>> but, as I mentioned before, the optimal timestep that gives the error
>>> bounds I need is just too small for my application, and the only other
>>> thing I can think to try is to throw a lot more cores at the problem and
>>> wait.
>>>
>>> Ellen
>>>
>>>
>>> On 10/15/19 4:14 PM, Smith, Barry F. wrote:
>>>>
>>>>  So you have a fixed "mesh" and fixed number of degrees of freedom for the 
>>>> entire time but the dependency on the function value at each particular 
>>>> point on the neighbor points changes over time?
>>>>
>>>>  For example, if you are doing upwinding and the direction changes when 
>>>> you used to use values from the right you now use values from the left?
>>>>
>>>>   Independent of the coloring, just changing the locations in the matrix

Re: [petsc-users] Changing nonzero structure and Jacobian coloring

[Ellen M. Price via petsc-users]

Thanks for the reply, Barry! Unfortunately, this is a particle code
(SPH, specifically), so the particle neighbors, which influence the
properties, change over time; the degrees of freedom is a constant, as
is the particle number. Any thoughts, given the new info? Or should I
stick with explicit integration? I've seen explicit used most commonly,
but, as I mentioned before, the optimal timestep that gives the error
bounds I need is just too small for my application, and the only other
thing I can think to try is to throw a lot more cores at the problem and
wait.

Ellen


On 10/15/19 4:14 PM, Smith, Barry F. wrote:
> 
>   So you have a fixed "mesh" and fixed number of degrees of freedom for the 
> entire time but the dependency on the function value at each particular point 
> on the neighbor points changes over time?
> 
>   For example, if you are doing upwinding and the direction changes when you 
> used to use values from the right you now use values from the left?
> 
>Independent of the coloring, just changing the locations in the matrix 
> where the entries are nonzero is expensive and painful. So what I would do is 
> build the initial Jacobian nonzero structure to contain all possible 
> connections, color that and then use that for the entire computation. At each 
> time step you will be working with some zero entries in the Jacobian but that 
> is ok, it is simpler and ultimately probably faster than trying to keep 
> changing the nonzero structure to "optimize" and only treat truly nonzero 
> values. 
> 
>For "stencil" type discretizations (finite difference, finite value) what 
> I describe should be fine. But if you problem is completely different (I 
> can't imagine how) and the Jacobian changes truly dramatically in structure 
> then what I suggest may not be appropriate but you'll need to tell us your 
> problem in that case so we can make suggestions.
> 
>Barry
> 
> 
> 
>> On Oct 15, 2019, at 2:56 PM, Ellen M. Price via petsc-users 
>>  wrote:
>>
>> Hi PETSc users!
>>
>> I have a problem that I am integrating with TS, and I think an implicit
>> method would let me take larger timesteps. The Jacobian is difficult to
>> compute, but, more importantly, the nonzero structure is changing with
>> time, so even if I use coloring and finite differences to compute the
>> actual values, I will have to update the coloring every time the
>> Jacobian recomputes.
>>
>> Has anyone done this before, and/or is there a correct way to tell TS to
>> re-compute the coloring of the matrix before SNES actually computes the
>> entries? Is this even a good idea, or is the coloring so expensive that
>> this is foolish (in general -- I know the answer depends on the problem,
>> but there may be a rule of thumb)?
>>
>> Thanks,
>> Ellen Price
> 

[petsc-users] Changing nonzero structure and Jacobian coloring

[Ellen M. Price via petsc-users]

Hi PETSc users!

I have a problem that I am integrating with TS, and I think an implicit
method would let me take larger timesteps. The Jacobian is difficult to
compute, but, more importantly, the nonzero structure is changing with
time, so even if I use coloring and finite differences to compute the
actual values, I will have to update the coloring every time the
Jacobian recomputes.

Has anyone done this before, and/or is there a correct way to tell TS to
re-compute the coloring of the matrix before SNES actually computes the
entries? Is this even a good idea, or is the coloring so expensive that
this is foolish (in general -- I know the answer depends on the problem,
but there may be a rule of thumb)?

Thanks,
Ellen Price

Re: [petsc-users] Using DMCOMPOSITE with TS

[Ellen M. Price via petsc-users]

Quick update: I found that changing TS_EXACTFINALTIME_INTERPOLATE to
TS_EXACTFINALTIME_MATCHSTEP fixes the problem; is there a reason I can't
use the interpolation when using a DMCOMPOSITE?

Ellen Price


On 2/20/19 1:42 PM, Ellen Price wrote:
> Hi fellow PETSc users!
> 
> I am attempting to use a DMCOMPOSITE alongside TS and have run into some
> trouble. I'm attaching a MWE demonstrating the problem. The goal is to
> combine a DMDA3d for spatial data and a DMREDUNDANT for non-spatial,
> time-dependent fields. In the attached example, this additional field is
> temperature.
> 
> The DMDA data is currently just temperature-dependent, and the
> temperature is supposed to increase linearly. Unfortunately, no matter
> how I integrate (explicitly, using Euler or RK, or implicitly, using
> backward Euler or BDF), I get the wrong final temperature. The
> integration proceeds for 100 timesteps and then stops (verified by
> -ts_monitor). At a heating rate of 1, and an initial temperature of 150,
> I should get a final temperature of 250 very easily. However, I get
> something closer to 200 (not exact, which makes me think this isn't a
> simple missing-a-factor-of-2 error).
> 
> I'm currently using TSSetDM with the composite DM to let PETSc compute
> the Jacobian. Having checked all the inputs and outputs as well as I
> can, this is the only step that seems like it may be causing a problem;
> yet even explicit integration doesn't work, and that shouldn't need the
> Jacobian at all. I'm at a loss.
> 
> Run using:
> ./mwe -implicit yes -ts_type beuler -ts_monitor OR ./mwe -implicit no
> -ts_type euler -ts_monitor
> 
> Thanks in advance for any help,
> Ellen Price

Re: [petsc-users] Increasing norm with finer mesh

[Ellen M. Price]

Maybe a stupid suggestion, but sometimes I forget to call the
*SetFromOptions function on my object, and then get confused when
changing the options has no effect. Just a thought from a fellow grad
student.

Ellen


On 10/16/2018 09:36 PM, Matthew Knepley wrote:
> On Tue, Oct 16, 2018 at 9:14 PM Weizhuo Wang  > wrote:
> 
> I just tried both, neither of them make a difference. I got exactly
> the same curve with either combination.
> 
> 
> I have a hard time believing you. If you make the residual tolerance
> much finer, your error will definitely change.
> I run tests every day that do exactly this. You can run them too, since
> they are just examples.
> 
>   Thanks,
> 
>      Matt
>  
> 
> Thanks!
> 
> Wang weizhuo
> 
> On Tue, Oct 16, 2018 at 8:06 PM Matthew Knepley  > wrote:
> 
> On Tue, Oct 16, 2018 at 7:26 PM Weizhuo Wang
> mailto:weizh...@illinois.edu>> wrote:
> 
> Hello again!
> 
> After some tweaking the code is giving right answers now.
> However it start to disagree with MATLAB results
> ('traditional' way using matrix inverse) when the grid is
> larger than 100*100. My PhD advisor and I suspects that the
> default dimension of the Krylov subspace is 100 in the test
> case we are running. If so, is there a way to increase the
> size of the subspace?
> 
> 
> 1) The default subspace size is 30, not 100. You can increase
> the subspace size using
> 
>        -ksp_gmres_restart n
> 
> 2) The problem is likely your tolerance. The default solver
> tolerance is 1e-5. You can change it using
> 
>        -ksp_rtol 1e-9
> 
>   Thanks,
> 
>      Matt
>  
> 
> 
> Disagrees.png
> 
> Thanks!
> 
> Wang Weizhuo
> 
> On Tue, Oct 9, 2018 at 2:50 AM Mark Adams  > wrote:
> 
> To reiterate what Matt is saying, you seem to have the
> exact solution on a 10x10 grid. That makes no sense
> unless the solution can be represented exactly by your
> FE space (eg, u(x,y) = x + y).
> 
> On Mon, Oct 8, 2018 at 9:33 PM Matthew Knepley
> mailto:knep...@gmail.com>> wrote:
> 
> On Mon, Oct 8, 2018 at 9:28 PM Weizhuo Wang
>  > wrote:
> 
> The code is attached in case anyone wants to
> take a look, I will try the high frequency
> scenario later.
> 
> 
> That is not the error. It is superconvergence at the
> vertices. The real solution is trigonometric, so your
> linear interpolants or whatever you use is not going
> to get the right value in between mesh points. You
> need to do a real integral over the whole interval
> to get the L_2 error.
> 
>   Thanks,
> 
>      Matt
>  
> 
> On Mon, Oct 8, 2018 at 7:58 PM Mark Adams
> mailto:mfad...@lbl.gov>> wrote:
> 
> 
> 
> On Mon, Oct 8, 2018 at 6:58 PM Weizhuo Wang
>  > wrote:
> 
> The first plot is the norm with the flag
> -pc_type lu with respect to number of
> grids in one axis (n), and the second
> plot is the norm without the flag
> -pc_type lu. 
> 
> 
> So you are using the default PC w/o LU. The
> default is ILU. This will reduce high
> frequency effectively but is not effective
> on the low frequency error. Don't expect
> your algebraic error reduction to be at the
> same scale as the residual reduction (what
> KSP measures). 
>  
> 
> 
> 
> -- 
> Wang Weizhuo
> 
> 
> 
> -- 
> 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 project function error

[Ellen M. Price]

Hi Matt,

I ran into a problem implementing this. My integration code relies on
knowing the points of each tetrahedron in the domain so I can turn an
integral over all space into a sum of integrals over cells. I can't find
an example of getting vertices for a given cell (I think TS ex11 comes
closest, looks like it gets faces?). Is there an elegant way to do this,
or an example I can follow?

Ellen


On 10/13/2018 04:26 AM, Matthew Knepley wrote:
> On Fri, Oct 12, 2018 at 6:00 PM Ellen M. Price
> mailto:ellen.pr...@cfa.harvard.edu>> wrote:
> 
> Hi Matt,
> 
> Thanks for the feedback! I'm always open to editorial comments. Do you
> have a suggestion for a better way? I'm basically just trying to do an
> integral over all space to get my boundary values for the FEM solve. I
> don't know if PETSc has a way to do something like this?
> 
> 
> So you are defining your BC by a boundary integral. I would just stick
> that integration code
> into a function and let that be called as your BC function. It will not
> be called "extra" times,
> just once for each dual vector, which in the case of P1 is every vertex.
> This way it will work
> for any element.
> 
>   Thanks,
> 
>      Matt
>  
> 
> Ellen
> 
> 
> On 10/12/2018 04:36 PM, Matthew Knepley wrote:
> > On Fri, Oct 12, 2018 at 3:51 PM Ellen M. Price
> > mailto:ellen.pr...@cfa.harvard.edu>
> <mailto:ellen.pr...@cfa.harvard.edu
> <mailto:ellen.pr...@cfa.harvard.edu>>> wrote:
> >
> >     Hi Matt,
> >
> >     Basically, I have a simple meshing script that generates a set
> of points
> >     for my FEM mesh and outputs a C file that compiles into my
> program.
> >     Along with those points is a static data field that I need to be
> >     associated with the correct data point. I don't have a
> function that
> >     evaluates at each point, it's just pre-generated data.
> >
> >
> > I cannot help but make editorial comments. Do not do this. It is
> > incredibly fragile, since
> > it depends on both a mesh and a discretization.
> >  
> >
> >     Since I called DMPlexDistribute() early on, I think my mesh
> points get
> >     re-ordered, right? So I need a way to re-order my data field
> to match
> >     that ordering. Is there a good way to do this?
> >
> >
> > However, if you want to do this, it is not hard. You do use
> > DMPlexDistributeField
> >
> 
> (https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DMPLEX/DMPlexDistributeField.html),
> > but you need
> > to be careful about the arguments. You want
> >
> >    DMPlexDistribute(dmOriginal, 0, &migrationSF, &dmNew);
> >    DMPlexDistributeField(dmOriginal, migrationSF, originalSection,
> > originalVec, newSection, newVec);
> >
> > Your Sections are easy to make, just put 1 dof on every vertex. The
> > originalVec is just a Vec of your values
> > (you can use VecCreateSeqWithArray() to wrap one around it), and the
> > newVec you can get with
> > DMCreateLocalVector(dmNew, &newVec) if you do DMSetSection(dmNew,
> > newSection).
> >
> > Does this make sense?
> >
> >   Thanks,
> >
> >      Matt
> >  
> >
> >     I tried using DMPlexDistributeField as a first attempt, but I
> get some
> >     kind of error no matter what I do (presumably I'm not giving
> the right
> >     input), and there aren't any examples in the docs.
> >
> >     Ellen
> >
> >
> >     On 10/12/2018 11:31 AM, Matthew Knepley wrote:
> >     > On Fri, Oct 12, 2018 at 11:12 AM Ellen M. Price
> >     >  <mailto:ellen.pr...@cfa.harvard.edu>
> <mailto:ellen.pr...@cfa.harvard.edu
> <mailto:ellen.pr...@cfa.harvard.edu>>
> >     <mailto:ellen.pr...@cfa.harvard.edu
> <mailto:ellen.pr...@cfa.harvard.edu>
> >     <mailto:ellen.pr...@cfa.harvard.edu
> <mailto:ellen.pr...@cfa.harvard.edu>>>> wrote:
> >     >
> >     >     So I found *that* problem -- turns out that using
> >     DMPlexCreateSection
> >     >     was messing things up, and it seems PETSc is handling
> that part
> >     >     internally?
> >     >
> >     >
> > 

Re: [petsc-users] DMPLEX project function error

[Ellen M. Price]

Hi Matt,

Basically, I have a simple meshing script that generates a set of points
for my FEM mesh and outputs a C file that compiles into my program.
Along with those points is a static data field that I need to be
associated with the correct data point. I don't have a function that
evaluates at each point, it's just pre-generated data.

Since I called DMPlexDistribute() early on, I think my mesh points get
re-ordered, right? So I need a way to re-order my data field to match
that ordering. Is there a good way to do this?

I tried using DMPlexDistributeField as a first attempt, but I get some
kind of error no matter what I do (presumably I'm not giving the right
input), and there aren't any examples in the docs.

Ellen


On 10/12/2018 11:31 AM, Matthew Knepley wrote:
> On Fri, Oct 12, 2018 at 11:12 AM Ellen M. Price
> mailto:ellen.pr...@cfa.harvard.edu>> wrote:
> 
> So I found *that* problem -- turns out that using DMPlexCreateSection
> was messing things up, and it seems PETSc is handling that part
> internally? 
> 
> 
> DMPlexSection creates the Section based on the discretization
> information in the DM (stored in a PetscDS).
> If no discretization info is available, then it defaults to P0 I think,
> which is what you got.
> 
> After discretization info is put into the DM, if I ask for the Section,
> the DM will create it automatically if it is missing.
>  
> 
> However, I've run into a related problem now. I need to
> project some field data onto the mesh; it's a static array of
> integration weights, and each one corresponds to a single point in the
> mesh. Based on this answer:
> 
> https://lists.mcs.anl.gov/pipermail/petsc-users/2015-December/027964.html
> 
> I think I need DMPlexDistributeField.
> 
> 
> I don't think that is what you want.  That is just sending data around.
> 
> You want to define an FEM field. We use DMProject*() for this. These
> take as input either
> a function of the coordinates (regular function) or another FEM field.
> 
> Can you tell me how your function is defined? I mean apart from the
> particular data you
> have at points.
> 
>   Thanks,
> 
>      Matt
>  
> 
> However, *that* function takes a
> section as an argument, and I had to take out the part where I create a
> section to get things working!  Calling DMGetDefaultSection or
> DMGetDefaultGlobalSection doesn't seem to help matters, either, because
> then I seem to get a null section?
> 
> So I either need to figure out how to set the section properly or how to
> get the section properly. As always, any help is appreciated.
> 
> Ellen
> 
> 
> On 10/11/2018 05:41 PM, Matthew Knepley wrote:
> > On Thu, Oct 11, 2018 at 4:39 PM Ellen M. Price
> > mailto:ellen.pr...@cfa.harvard.edu>
> <mailto:ellen.pr...@cfa.harvard.edu
> <mailto:ellen.pr...@cfa.harvard.edu>>> wrote:
> >
> >     I was working with a DMPLEX and FEM, following SNES example
> 12. I get
> >     the following error when I call DMProjectFunction, but I don't
> know what
> >     it means. Can anyone explain where I might have gone wrong, or
> at least
> >     what this error is telling me? I think the point closure size is
> >     correct, since my mesh is 3d simplex,
> >
> >
> > Yes, if you have 3D SIMPLEX mesh and are using P1 elements, then you
> > would have
> > 4 dofs in the closure of a cell. The dual space dimension is the
> number
> > of dual space
> > basis vectors assigned to points in the closure. Since it is 1, it
> looks
> > like you have a P0
> > dual space. I assume you changed something in ex12?
> >
> >   Thanks,
> >
> >      Matt
> >  
> >
> >     but what is the dual space
> >     dimension, and where might I have set it incorrectly?
> >
> >     [0]PETSC ERROR: Nonconforming object sizes
> >     [0]PETSC ERROR: The section point closure size 4 != dual space
> >     dimension 1
> >     [0]PETSC ERROR: See
> http://www.mcs.anl.gov/petsc/documentation/faq.html
> >     for trouble shooting.
> >     [0]PETSC ERROR: Petsc Release Version 3.9.2, May, 20, 2018
> >     ...
> >     [0]PETSC ERROR: Configure options
> >     --prefix=/home/eprice/software/petsc-opt --with-hdf5=1
> >     --with-hdf5-dir=/home/eprice/software/hdf5-parallel --with-mpe=1
> >     --with-mpe-dir=/home/eprice/software/mpe --with-debugging=0
> >  

Re: [petsc-users] DMPLEX project function error

[Ellen M. Price]

So I found *that* problem -- turns out that using DMPlexCreateSection
was messing things up, and it seems PETSc is handling that part
internally? However, I've run into a related problem now. I need to
project some field data onto the mesh; it's a static array of
integration weights, and each one corresponds to a single point in the
mesh. Based on this answer:

https://lists.mcs.anl.gov/pipermail/petsc-users/2015-December/027964.html

I think I need DMPlexDistributeField. However, *that* function takes a
section as an argument, and I had to take out the part where I create a
section to get things working!  Calling DMGetDefaultSection or
DMGetDefaultGlobalSection doesn't seem to help matters, either, because
then I seem to get a null section?

So I either need to figure out how to set the section properly or how to
get the section properly. As always, any help is appreciated.

Ellen


On 10/11/2018 05:41 PM, Matthew Knepley wrote:
> On Thu, Oct 11, 2018 at 4:39 PM Ellen M. Price
> mailto:ellen.pr...@cfa.harvard.edu>> wrote:
> 
> I was working with a DMPLEX and FEM, following SNES example 12. I get
> the following error when I call DMProjectFunction, but I don't know what
> it means. Can anyone explain where I might have gone wrong, or at least
> what this error is telling me? I think the point closure size is
> correct, since my mesh is 3d simplex,
> 
> 
> Yes, if you have 3D SIMPLEX mesh and are using P1 elements, then you
> would have
> 4 dofs in the closure of a cell. The dual space dimension is the number
> of dual space
> basis vectors assigned to points in the closure. Since it is 1, it looks
> like you have a P0
> dual space. I assume you changed something in ex12?
> 
>   Thanks,
> 
>      Matt
>  
> 
> but what is the dual space
> dimension, and where might I have set it incorrectly?
> 
> [0]PETSC ERROR: Nonconforming object sizes
> [0]PETSC ERROR: The section point closure size 4 != dual space
> dimension 1
> [0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html
> for trouble shooting.
> [0]PETSC ERROR: Petsc Release Version 3.9.2, May, 20, 2018
> ...
> [0]PETSC ERROR: Configure options
> --prefix=/home/eprice/software/petsc-opt --with-hdf5=1
> --with-hdf5-dir=/home/eprice/software/hdf5-parallel --with-mpe=1
> --with-mpe-dir=/home/eprice/software/mpe --with-debugging=0
> LDFLAGS="-pthread -lz" COPTFLAGS="-O3 -march=native -mtune=native"
> CXXOPTFLAGS="-O3 -march=native -mtune=native" FOPTFLAGS="-O3
> -march=native -mtune=native" --with-mpi=1
> --with-mpi-dir=/home/eprice/software/mpich --with-mumps=1
> --with-mumps-dir=/home/eprice/software/mumps --with-parmetis=1
> --with-parmetis-dir=/home/eprice/software/parmetis --with-metis=1
> --with-metis-dir=/home/eprice/software/parmetis --with-ptscotch=1
> --with-ptscotch-dir=/home/eprice/software/scotch --with-scalapack=1
> --with-scalapack-dir=/home/eprice/software/scalapack
> [0]PETSC ERROR: #1 DMProjectLocal_Generic_Plex() line 347 in
> /h/sabriel0/src/petsc-3.9.2/src/dm/impls/plex/plexproject.c
> [0]PETSC ERROR: #2 DMProjectFunctionLocal_Plex() line 428 in
> /h/sabriel0/src/petsc-3.9.2/src/dm/impls/plex/plexproject.c
> [0]PETSC ERROR: #3 DMProjectFunctionLocal() line 6265 in
> /h/sabriel0/src/petsc-3.9.2/src/dm/interface/dm.c
> [0]PETSC ERROR: #4 DMProjectFunction() line 6250 in
> /h/sabriel0/src/petsc-3.9.2/src/dm/interface/dm.c
> ...
> 
> (I know this is an optimized PETSc build, but I get the same error from
> my debug build, it's just much slower.)
> 
> Ellen
> 
> 
> 
> -- 
> What most experimenters take for granted before they begin their
> experiments is infinitely more interesting than any results to which
> their experiments lead.
> -- Norbert Wiener
> 
> https://www.cse.buffalo.edu/~knepley/
> <http://www.cse.buffalo.edu/%7Eknepley/>

[petsc-users] DMPLEX project function error

[Ellen M. Price]

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

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

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

Ellen

Re: [petsc-users] Implementing FEM with PLEX/DS/SNES

[Ellen M. Price]

Hi Matt,

Thanks for the speedy reply. I do have a few more questions/elaborations
on my previous correspondence.

As for the source term (and the boundary condition, actually), both are
only *known* at mesh points, but are assumed continuous everywhere. How
does this change things? Basically, I want to assume a linear
interpolation (I'm using linear elements) of the source term over each
tetrahedral cell and then the integral should be some linear combination
of mesh point values, with coefficients given by the inner product of
basis functions? The boundary term is actually given by a very nasty
integral equation, which is probably the most expensive part of the
simulation. So I only want to evaluate it at mesh points if possible.

You are correct, I forgot the \cdot n in that equation. Oops!

I need to think more about the Jacobian and residual before I can really
ask an intelligent question about them. I suppose that in practice, I
need to figure out exactly which Jacobian terms I need to include, and
which can be NULL. Some guidance on that would be immensely helpful.

If it helps, the equation I'm solving is just the Newtonian
gravitational potential equation (on a nontrivial mesh and with a
nontrivial, non-analytic mass density function).

Ellen


On 10/08/2018 06:28 PM, Matthew Knepley wrote:
> On Mon, Oct 8, 2018 at 5:50 PM Ellen M. Price
> mailto:ellen.pr...@cfa.harvard.edu>> wrote:
> 
> Can anyone shed some light on SNES example 12?
> 
> 
> https://www.mcs.anl.gov/petsc/petsc-current/src/snes/examples/tutorials/ex12.c.html
> 
> I am solving a Poisson problem, but with a source term defined on mesh
> points, so I don't have an analytic form.
> 
> 
> Can you explain more about this source term? It sounds like a bunch of
> delta functions. That
> would still work in this framework, but the convergence rate for a rhs
> with these singularities
> is reduced (this is a generic feature of FEM).
>  
> 
> The part I am confused about
> is the role of the functions f0, f1, and g3 that go into the DS.
> 
> Suppose I have the Poisson problem nabla^2 u = f. Then the weak form is,
> for test function phi,
> 
> \int_{\partial \Omega} \phi \nabla u - \int_\Omega \nabla u \cdot \nabla
> \phi = \int_\Omega \phi f
> 
> 
> I think you need \int_{\partial \Omega} \phi \nabla u \cdot n instead.
>  
> 
> Based on the documentation of PetscDSSetResidual(), I'm a little
> confused about f0 and f1:
> 
> \int_\Omega \phi f_0(u, u_t, \nabla u, x, t) + \nabla\phi \cdot {\vec
> f}_1(u, u_t, \nabla u, x, t)
> 
> So f0 is, for the Poisson problem, zero, and f1 is just -grad u?
> 
> 
> No. f1 is indeed -grad u, but f0 is -f.
> 
> Note that ex12 solves -nabla^2 u = f since that operator is positive
> definite.
>  
> 
> Do we
> need to do anything with the surface term, or is that handled
> internally?
> 
> 
> The surface term can be included using
> 
>   
> https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DT/PetscDSSetBdResidual.html
> 
> and we use that in ex12 for inhomogeneous Neumann conditions.
>  
> 
> Because neither of these terms are a surface integral term,
> both are volume integrals.
> 
> I can't figure out the purpose of g0, g1, g2, g3 in
> PetscDSSetJacobian(), either. Is g supposed to be the source term?
> 
> 
> No, it is the coefficient of the different combinations of basis
> functions and their gradients.
> For instance, the g3 term (grad phi and grad psi), g3 plays the same
> role as the elasticity
> tensor C in mechanics.
>  
> 
> Is
> psi a second basis function?
> 
> 
> Yes. If the method is Galerkin, it comes from the same space as phi.
>  
> 
> What even *is* the Jacobian in a FEM
> problem?
> 
> 
> The Frechet derivative of the residual.
>  
> 
> I've only seen them formulated as linear problems in tutorials,
> so I'm not even clear why we need a nonlinear solver or Jacobian...
> 
> 
> The Frechet derivative of a linear operator is itself. ex12 also
> incorporates
> nonlinear variants of Poisson.
>  
> 
> Any help is appreciated. I like the idea of a general framework for FEM,
> as I think it will reduce human error on my part, but I also want to
> understand what I'm doing.
> 
> 
> If you have any more questions, feel free to keep mailing.
> 
>   Thanks,
> 
>      Matt
>  
> 
> Ellen Price
> 
> 
> 
> -- 
> 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/
> <http://www.cse.buffalo.edu/%7Eknepley/>

[petsc-users] Implementing FEM with PLEX/DS/SNES

[Ellen M. Price]

Can anyone shed some light on SNES example 12?

https://www.mcs.anl.gov/petsc/petsc-current/src/snes/examples/tutorials/ex12.c.html

I am solving a Poisson problem, but with a source term defined on mesh
points, so I don't have an analytic form. The part I am confused about
is the role of the functions f0, f1, and g3 that go into the DS.

Suppose I have the Poisson problem nabla^2 u = f. Then the weak form is,
for test function phi,

\int_{\partial \Omega} \phi \nabla u - \int_\Omega \nabla u \cdot \nabla
\phi = \int_\Omega \phi f

Based on the documentation of PetscDSSetResidual(), I'm a little
confused about f0 and f1:

\int_\Omega \phi f_0(u, u_t, \nabla u, x, t) + \nabla\phi \cdot {\vec
f}_1(u, u_t, \nabla u, x, t)

So f0 is, for the Poisson problem, zero, and f1 is just -grad u? Do we
need to do anything with the surface term, or is that handled
internally? Because neither of these terms are a surface integral term,
both are volume integrals.

I can't figure out the purpose of g0, g1, g2, g3 in
PetscDSSetJacobian(), either. Is g supposed to be the source term? Is
psi a second basis function? What even *is* the Jacobian in a FEM
problem? I've only seen them formulated as linear problems in tutorials,
so I'm not even clear why we need a nonlinear solver or Jacobian...

Any help is appreciated. I like the idea of a general framework for FEM,
as I think it will reduce human error on my part, but I also want to
understand what I'm doing.

Ellen Price

[petsc-users] Question about HDF5 viewer groups

[Ellen M. Price]

Hi there,

I'm running a complex PETSc program that outputs a LOT of data to an
HDF5 file. As such, I'd like to separate my output into groups to make
it easier to traverse in the analysis phase.

I'm either confused about the nomenclature or the usage of
PetscViewerHDF5PushGroup/PetscViewerHDF5PopGroup. From a quick skim of
the source, it looks like the PushGroup function builds up a linked list
of group names? Also, just the name suggests that there should be a
"stack" of group names stored somewhere.

So I would expect, after calling PetscViewerHDF5PushGroup twice, say,

PetscViewerHDF5PushGroup(h5viewer, "group1");
PetscViewerHDF5PushGroup(h5viewer, "group2");
VecView(x, h5viewer);

that I would get my vector output in /group1/group2/x. This is obviously
not what's actually happening. If I run
http://www.mcs.anl.gov/petsc/petsc-current/src/vec/vec/examples/tutorials/ex19.c,
for example, I get, from h5dump:

HDF5 "ex19.h5" {
FILE_CONTENTS {
 group  /
 dataset/TestVec
 group  /testBlockSize
 dataset/testBlockSize/TestVec2
 group  /testTimestep
 dataset/testTimestep/TestVec2
 }
}

even though multiple groups are "pushed" in sequence.

My question is, have I misunderstood the connotation of pushing/popping
groups? How can I properly create a group and then a subgroup? Every
attempt at doing this:

PetscViewerHDF5PushGroup(h5viewer, "/group1");
PetscViewerHDF5PushGroup(h5viewer, "/group1/group2");
VecView(x, h5viewer);

leads to a long string of HDF5 library errors, like:

HDF5-DIAG: Error detected in HDF5 (1.8.16) MPI-process 0:
  #000: ../../../src/H5L.c line 824 in H5Lexists(): unable to get link info
major: Symbol table
minor: Object not found
  #001: ../../../src/H5L.c line 2765 in H5L_exists(): path doesn't exist
major: Symbol table
minor: Object already exists
  #002: ../../../src/H5Gtraverse.c line 861 in H5G_traverse(): internal
path traversal failed
major: Symbol table
minor: Object not found
  #003: ../../../src/H5Gtraverse.c line 755 in H5G_traverse_real():
component not found
major: Symbol table
minor: Object not found
HDF5-DIAG: Error detected in HDF5 (1.8.16) MPI-process 0:
  #000: ../../../src/H5G.c line 314 in H5Gcreate2(): unable to create group
major: Symbol table
minor: Unable to initialize object
  #001: ../../../src/H5Gint.c line 194 in H5G__create_named(): unable to
create and link to group
major: Symbol table
minor: Unable to initialize object
  #002: ../../../src/H5L.c line 1638 in H5L_link_object(): unable to
create new link to object
major: Links
minor: Unable to initialize object
  #003: ../../../src/H5L.c line 1882 in H5L_create_real(): can't insert link
major: Symbol table
minor: Unable to insert object
  #004: ../../../src/H5Gtraverse.c line 861 in H5G_traverse(): internal
path traversal failed
major: Symbol table
minor: Object not found
  #005: ../../../src/H5Gtraverse.c line 755 in H5G_traverse_real():
component not found
major: Symbol table
minor: Object not found

Thanks in advance,
Ellen Price

[petsc-users] Debugging SNES when "everything" looks okay

[Ellen M. Price]

Hi PETSc users,

I'm having trouble applying SNES to a new problem I'm working on. I'll
try to be as complete as possible but can't post the full code because
it's ongoing research and is pretty long anyway.

The nonlinear problem arises from trying to solve a set of two coupled
ODEs using a Galerkin method. I am using Mathematica to generate the
objective function to solve and the Jacobian, so I *think* I can rule
out human error on that front.

There are four things I can easily change:

- number of DMDA grid points N (I've tried 100 and 1000)
- preconditioner (I've tried LU and SVD, LU appears to work better, and
SVD is too slow for N = 1000)
- linear solver (haven't played with this much, but sometimes smaller
tolerances work better)
- nonlinear solver (what I'm having trouble with)

Trying different solvers should, as far as I'm aware, give comparable
answers, but that's not the case here. NEWTONTR works best of the ones
I've tried, but I'm suspicious that the function value barely decreases
before SNES "converges" -- and none of the options I've tried changing
seem to affect this, as it just finds another reason to converge without
making any real progress. For example:

  0 SNES Function norm 7.197788479418e+00
0 KSP Residual norm 1.721996766619e+01
1 KSP Residual norm 5.186021549059e-14
  Linear solve converged due to CONVERGED_STEP_LENGTH iterations 1
  1 SNES Function norm 7.19674987e+00
0 KSP Residual norm 3.296112185897e+01
1 KSP Residual norm 2.713415432045e-13
  Linear solve converged due to CONVERGED_STEP_LENGTH iterations 1
.
 50 SNES Function norm 7.197376803542e+00
0 KSP Residual norm 6.222518656302e+02
1 KSP Residual norm 9.630659996504e-12
  Linear solve converged due to CONVERGED_STEP_LENGTH iterations 1
 51 SNES Function norm 7.197376803542e+00
Nonlinear solve converged due to CONVERGED_SNORM_RELATIVE iterations 50

I've tried everything I can think of and the FAQ suggestions, including
non-dimensionalizing the problem; I observe the same behavior either
way. The particularly strange thing I can't understand is why some of
the SNES methods fail outright, after just one iteration (NCG, NGMRES,
and others) with DIVERGED_DTOL. Unless I've misunderstood, it seems
like, for the most part, I should be able to substitute in one method
for another, possibly adjusting a few parameters along the way.
Furthermore, the default method, NEWTONLS, diverges with
DIVERGED_LINE_SEARCH, which I'm not sure how to debug.

So the only viable method is NEWTONTR, and that doesn't appear to
"really" converge. Any suggestions for further things to try are
appreciated. My current options are:

-pc_type lu -snes_monitor -snes_converged_reason -ksp_converged_reason
-snes_max_it 1 -ksp_monitor -snes_type newtonls

Thanks in advance,
Ellen Price

[petsc-users] Errors encountered when linking against a successful build including MUMPS

[Ellen M. Price]

I just successfully (re)-compiled petsc-3.9.2 to use my MUMPS install.
The PETSc build process completed with no errors and I installed my new
shared libs. However, when I tried linking against the new libraries, I
get the following errors:

/home/eprice/software/petsc-opt/lib/libpetsc.so: undefined reference to
`MatMumpsSetIcntl'
/home/eprice/software/petsc-opt/lib/libpetsc.so: undefined reference to
`MatPartitioningParmetisSetRepartition'
/home/eprice/software/petsc-opt/lib/libpetsc.so: undefined reference to
`MatSolverTypeRegister_MUMPS'
/home/eprice/software/petsc-opt/lib/libpetsc.so: undefined reference to
`MatPartitioningCreate_Parmetis'
/home/eprice/software/petsc-opt/lib/libpetsc.so: undefined reference to
`MatPartitioningCreate_PTScotch'
collect2: error: ld returned 1 exit status

This is confusing, because I can't see any reason these functions
shouldn't be included in the final build. My configure command follows:

./configure --prefix=/home/eprice/software/petsc-opt --with-hdf5=1
--with-hdf5-dir=/home/eprice/software/hdf5-parallel --with-mpe=1
--with-mpe-dir=/home/eprice/software/mpe --with-debugging=0
LDFLAGS='-pthread -lz' COPTFLAGS='-O3 -march=native -mtune=native'
CXXOPTFLAGS='-O3 -march=native -mtune=native' FOPTFLAGS='-O3
-march=native -mtune=native' --with-mpi=1
--with-mpi-dir=/home/eprice/software/mpich --with-mumps=1
--with-mumps-dir=/home/eprice/software/mumps --with-parmetis=1
--with-parmetis-dir=/home/eprice/software/parmetis --with-metis=1
--with-metis-dir=/home/eprice/software/parmetis --with-ptscotch=1
--with-ptscotch-dir=/home/eprice/software/scotch --with-scalapack=1
--with-scalapack-dir=/home/eprice/software/scalapack

Any help is greatly appreciated. I tried picking through the build
system to find a problem, but it's difficult to understand without docs.

Ellen

Re: [petsc-users] What does DMPlexDistribute actually do?

[Ellen M. Price]

Here is the full minimal (non)-working example and the problematic
output. I'm using PETSc 3.8.0.


#include 

#include 


int main(int argc, char *argv[])
{
PetscInt refine = 0;
PetscErrorCode ierr;
DM dmplex, dmplex_dist;


/* initialize */
PetscInitialize(&argc, &argv, NULL, NULL);


/* get the refinement level and create the mesh */
PetscOptionsGetInt(NULL, NULL, "-refine", &refine, NULL);
ierr = DMPlexCreateHexCylinderMesh(PETSC_COMM_WORLD, refine,
DM_BOUNDARY_NONE, &dmplex); CHKERRQ(ierr);
DMView(dmplex, PETSC_VIEWER_STDOUT_WORLD);
ierr = DMPlexDistribute(dmplex, 0,
NULL, &dmplex_dist); CHKERRQ(ierr);
if (dmplex_dist)
{
DMDestroy(&dmplex);
dmplex = dmplex_dist;
}
DMView(dmplex, PETSC_VIEWER_STDOUT_WORLD);
DMDestroy(&dmplex);

/* finalize */

PetscFinalize();

return 0;
}


$ mpirun -n 2 bin/mwe -refine 2

DM Object: 2 MPI processes
  type: plex
DM_0x1bde850_0 in 3 dimensions:
  0-cells: 445 445
  1-cells: 1196 1196
  2-cells: 1072 1072
  3-cells: 320 320
Labels:
  depth: 4 strata with value/size (0 (445), 1 (1196), 2 (1072), 3 (320))
DM Object: Parallel Mesh 2 MPI processes
  type: plex
Parallel Mesh in 3 dimensions:
  0-cells: 445 445
  1-cells: 1196 1196
  2-cells: 1072 1072
  3-cells: 320 320
Labels:
  depth: 4 strata with value/size (0 (445), 1 (1196), 2 (1072), 3 (320))

$ mpirun -n 4 bin/mwe -refine 2

DM Object: 4 MPI processes
  type: plex
DM_0x20c7790_0 in 3 dimensions:
  0-cells: 445 445 445 445
  1-cells: 1196 1196 1196 1196
  2-cells: 1072 1072 1072 1072
  3-cells: 320 320 320 320
Labels:
  depth: 4 strata with value/size (0 (445), 1 (1196), 2 (1072), 3 (320))
DM Object: Parallel Mesh 4 MPI processes
  type: plex
Parallel Mesh in 3 dimensions:
  0-cells: 445 445 445 445
  1-cells: 1196 1196 1196 1196
  2-cells: 1072 1072 1072 1072
  3-cells: 320 320 320 320
Labels:
  depth: 4 strata with value/size (0 (445), 1 (1196), 2 (1072), 3 (320))


Ellen Price


On 03/25/2018 07:29 PM, Matthew Knepley wrote:
> On Sun, Mar 25, 2018 at 6:32 PM, Ellen M. Price
> mailto:ellen.pr...@cfa.harvard.edu>> wrote:
> 
> I am trying to understand some unusual behavior (at least, given my
> current understanding) in DMPlexDistribute. I have created a hex mesh
> and distributed it using the following snippet:
> 
> /* "-refine" is an option set at runtime */
> PetscOptionsGetInt(NULL, NULL, "-refine", &refine, NULL);
> ierr = DMPlexCreateHexCylinderMesh(PETSC_COMM_WORLD, refine,
>     DM_BOUNDARY_NONE, &dmplex); CHKERRQ(ierr);
> DMView(dmplex, PETSC_VIEWER_STDOUT_WORLD);
> 
> /* this is from the examples */
> ierr = DMPlexDistribute(dmplex, 0, NULL, &dmplex_dist); CHKERRQ(ierr);
> 
> if (dmplex_dist)
> 
> {
> 
>     DMDestroy(&dmplex);
> 
>     dmplex = dmplex_dist;
> 
> }
> 
> DMView(dmplex, PETSC_VIEWER_STDOUT_WORLD);
> 
> So I do a view operation before and after the distribute call to see how
> the DM is structured. I do not understand what happens next:
> 
> $ mpirun -n 4 ./myprogram -refine 2
> 
> DM Object: 4 MPI processes
>   type: plex
> DM_0x1f24d50_0 in 3 dimensions:
>   0-cells: 445 445 445 445
>   1-cells: 1196 1196 1196 1196
>   2-cells: 1072 1072 1072 1072
>   3-cells: 320 320 320 320
> Labels:
>   depth: 4 strata with value/size (0 (445), 1 (1196), 2 (1072), 3 (320))
> DM Object: Parallel Mesh 4 MPI processes
>   type: plex
> Parallel Mesh in 3 dimensions:
>   0-cells: 445 445 445 445
>   1-cells: 1196 1196 1196 1196
>   2-cells: 1072 1072 1072 1072
>   3-cells: 320 320 320 320
> Labels:
>   depth: 4 strata with value/size (0 (445), 1 (1196), 2 (1072), 3 (320))
> 
> No matter what I choose for the number of processors, every processor
> has a copy of all 320 cells (at this refinement level). Similar behavior
> at other refinement levels. The only thing that is changed is the name
> of the DM, to "Parallel Mesh". This is not what I would have expected
> given the description of DMPlexDistribute in the manual; I thought the
> cells would be split up between all available processors.
> 
> I am also viewing this mesh in VTK and have noticed that the file size
> of the output scales with the number of processors, as if it is really
> writing each "copy" of the mesh and data stored in it to one big file.
> Again, not what I expected.
> 
> Can someone clear up what is going on?
> 
> 
> It definitely does not work that way. Can you send the whole code you
> are running?
> Its hard to tell where the pro

[petsc-users] What does DMPlexDistribute actually do?

[Ellen M. Price]

I am trying to understand some unusual behavior (at least, given my
current understanding) in DMPlexDistribute. I have created a hex mesh
and distributed it using the following snippet:

/* "-refine" is an option set at runtime */
PetscOptionsGetInt(NULL, NULL, "-refine", &refine, NULL);
ierr = DMPlexCreateHexCylinderMesh(PETSC_COMM_WORLD, refine,
DM_BOUNDARY_NONE, &dmplex); CHKERRQ(ierr);
DMView(dmplex, PETSC_VIEWER_STDOUT_WORLD);

/* this is from the examples */
ierr = DMPlexDistribute(dmplex, 0, NULL, &dmplex_dist); CHKERRQ(ierr);

if (dmplex_dist)

{

DMDestroy(&dmplex);

dmplex = dmplex_dist;

}

DMView(dmplex, PETSC_VIEWER_STDOUT_WORLD);

So I do a view operation before and after the distribute call to see how
the DM is structured. I do not understand what happens next:

$ mpirun -n 4 ./myprogram -refine 2

DM Object: 4 MPI processes
  type: plex
DM_0x1f24d50_0 in 3 dimensions:
  0-cells: 445 445 445 445
  1-cells: 1196 1196 1196 1196
  2-cells: 1072 1072 1072 1072
  3-cells: 320 320 320 320
Labels:
  depth: 4 strata with value/size (0 (445), 1 (1196), 2 (1072), 3 (320))
DM Object: Parallel Mesh 4 MPI processes
  type: plex
Parallel Mesh in 3 dimensions:
  0-cells: 445 445 445 445
  1-cells: 1196 1196 1196 1196
  2-cells: 1072 1072 1072 1072
  3-cells: 320 320 320 320
Labels:
  depth: 4 strata with value/size (0 (445), 1 (1196), 2 (1072), 3 (320))

No matter what I choose for the number of processors, every processor
has a copy of all 320 cells (at this refinement level). Similar behavior
at other refinement levels. The only thing that is changed is the name
of the DM, to "Parallel Mesh". This is not what I would have expected
given the description of DMPlexDistribute in the manual; I thought the
cells would be split up between all available processors.

I am also viewing this mesh in VTK and have noticed that the file size
of the output scales with the number of processors, as if it is really
writing each "copy" of the mesh and data stored in it to one big file.
Again, not what I expected.

Can someone clear up what is going on?

Ellen Price