Re: [petsc-users] MatCompositeMerge + MatCreateRedundantMatrix

[Marius Buerkle via petsc-users]

ok, I think I understand now. I will give it a try and if there is some trouble comeback to you. thanks.

 

marius



 



On Tue, Feb 19, 2019 at 8:42 PM Marius Buerkle  wrote:




ok, so it seems there is no straight forward way to transfer data between PETSc matrices on different subcomms. Probably doing it by "hand" extracting the matricies on the subcomms create a MPI_INTERCOMM transfering the data to PETSC_COMM_WORLD and assembling them in a new PETSc matrix would be possible, right?



 

That sounds too complicated. Why not just reverse MPICreateSubMatricesMPI()? Meaning make it collective on the whole big communicator, so that you can swap out all the subcommunicator for the aggregation call, just like we do in that function.

Then its really just a matter of reversing the communication call.

 

   Matt 






 



On Tue, Feb 19, 2019 at 7:12 PM Marius Buerkle  wrote:





I see. This would work if the matrices are on different subcommumicators ? Is it possible to add this functionality ?




 

Hmm, no. That is specialized to serial matrices. You need the inverse of MatCreateSubMatricesMPI().

 

  Thanks,

 

     Matt

  




marius

 

 


You basically need the inverse of MatCreateSubmatrices(). I do not think we have that right now, but it could probably be done without too much trouble by looking at that code.
 

  Thanks,

 

     Matt

 


On Tue, Feb 19, 2019 at 6:15 AM Marius Buerkle via petsc-users  wrote:




Hi !

 

Is there some way to combine MatCompositeMerge with MatCreateRedundantMatrix? I basically want to create copies of a matrix from PETSC_COMM_WORLD to subcommunicators, do some work on each subcommunicator and than gather the results back to PETSC_COMM_WORLD, namely  I want to sum the  the invidual matrices from the subcommunicatos component wise and get the resulting matrix on PETSC_COMM_WORLD. Is this somehow possible without going through all the hassle of using MPI directly? 

 

marius




 

 
--







What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.
-- Norbert Wiener

 

https://www.cse.buffalo.edu/~knepley/













 

 
--







What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.
-- Norbert Wiener

 

https://www.cse.buffalo.edu/~knepley/














 

 
--







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] Including constrained dofs in a matrix/vector created with DMCreateXXXX()

[Jordan Wagner via petsc-users]

On 2/19/19 7:40 AM, Matthew Knepley wrote:
On Tue, Feb 19, 2019 at 1:13 AM Jordan Wagner via petsc-users 
mailto:petsc-users@mcs.anl.gov>> wrote:


Hi,

Over the past few months, I have implemented dmplex/section in my
preexisting finite element code. We already had our own
implementations for things like the finite element space, boundary
conditions, projections, etc, so I am mostly using dmplex/section
at its lowest level; no DT/DS stuff for the moment. Everything is
working fine at this level, and I am overall very happy with how
it ties into the code that we already have.

I am currently doing my own preallocation for the system Jacobian
but am wanting to use DMCreateMatrix() instead, as it seems
pointless to keep using my own crappy function when this one is
just sitting there waiting. However, I dislike how I am currently
implementing this and was hoping to get some suggestions or some
advice. My problem is handling nicely the Dirichlet dofs. My code
takes the approach of first assembling the entire Jacobian matrix
and load vector as if no constraints are imposed. We then have a
function that applies the essential conditions after assembly and
extracts a submatrix and subvector similar to the way ex49


is doing.

Since DMCreateMatrix() automatically partitions out the
constraints specified in the section, I find myself having to
create two nearly identical sections: one that has constraints and
one that does not, and then I clone the DM and set the default
section of each respectively. I can then use DMCreateMatrix() on
the one dm/section with no constraints to preallocate the larger
matrix that I am assembling. Then I use the dm/section with
constraints to do all of my boundary condition and projection
operations on the previously created matrix and vector.
Essentially, I am using an unconstrained section only for
allocation and a constrained section for most everything else.

This process actually works fine for me, but its pretty ugly, and
I am sure there is a better way. I am wondering if I am missing
something that could keep me from having to go through this
process of creating two sections that differ only in the
constraints. It seems to me if there were an option for
DMCreateMatrix()/DMCreateVector() to either include or not include
the dofs, that would completely solve my problem. That way, I can
use the same section for both creation and allocation.

So, the main question I have is just whether there is a function
or option that I am missing like DMCreate() but with the
ability to retain the constrained rows and columns rather than
compressing them out.

I do not completely understand the above, but here is how Section 
works. A _local_ Section contains all dofs and is intended for 
assembly using local Vecs. To communicate with the solver, you make a 
_global_ Section using PetscSectionCreateGlobalSection(). This has an 
option to include or exclude constraints. For example, I make a global 
Section with constraints when I output for visualization. Maybe you 
can use that to do what you want.


Thanks for the prompt reply and sorry for the unclear initial question. 
I think this helps in my understanding of how this tool was intended to 
be used.


What I do not understand is why you do not just use a global Section 
without constraints to assemble your matrix, rather than eliminating 
the rows/columns later. The global Section gives back negative indices 
for unknown that are constrained or not owned. These are ignored by 
MatSetValues() and VecSetValues(). This is how the default Plex 
assembly makes operators without constrained variables.


I have rarely used a global section directly because I have never been 
able to  make sense of exactly what it was giving me. For instance as 
you mentioned (and I have seen written elsewhere), the global section 
should give negative indices for unowned or constrained dofs. This seems 
to work fine for me with the unowned points---I clearly get a negative 
size and offset from the section on these points. But my problem is for 
points that have constraints. I don't see where I would get negative 
indices from. The sizes and offsets for these points are not negative.  
Of course, from the constraint size and constraint indices on these 
points, I can manually find out which should be negative and go from 
there. However, that leads me to the next problem I am having with 
global section...


For some reason, when I am creating the global section from the local 
one, my constrained component indices are getting messed up. The sizes 
are right but the order seems to be random and/or nonsensical. For 
instance, here is a snippet from my code being run in serial


Local section:

  (17615) 

[petsc-users] Using DMCOMPOSITE with TS

[Ellen Price via petsc-users]

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

#define ZERO(0.)

#define ONE (1.)

#define AUXTEMP (0)

#define TEMP0   (150.)

#define HEATING (1.)

#define NU  (0.1)

#define TFINAL  (100.)

#define DELTA   (1.)

#define unused(x)   ((void) (x))

PetscErrorCode rhsfunc_petsc(TS ts, PetscReal t, Vec y, Vec ydot, void *ptr);
PetscErrorCode ifunc_petsc(TS ts, PetscReal t, Vec y, Vec ydot, Vec f, void *ptr);
void func(double t, Vec y, Vec aux, Vec ydot, Vec auxdot, DM da);

PetscErrorCode rhsfunc_petsc(TS ts, PetscReal t, Vec y, Vec ydot, void *ptr)
{
DM pack, da;
Vec yda, yred, ydotda, ydotred;

unused(ptr);

TSGetDM(ts, );

DMCompositeGetEntries(pack, , NULL);

DMCompositeGetLocalVectors(pack, , );
DMCompositeGetLocalVectors(pack, , );
DMCompositeScatter(pack, y, yda, yred);

VecSet(ydotda, ZERO);
VecSet(ydotred, ZERO);

func(t, yda, yred, ydotda, ydotred, da);

DMCompositeGather(pack, INSERT_VALUES, ydot, ydotda, ydotred);
DMCompositeRestoreLocalVectors(pack, , );
DMCompositeRestoreLocalVectors(pack, , );

VecView(ydot, PETSC_VIEWER_STDOUT_WORLD);

return 0;
}

PetscErrorCode ifunc_petsc(TS ts, PetscReal t, Vec y, Vec ydot, Vec f, void *ptr)
{
DM pack, da;
Vec yda, yred, ydotda, ydotred;

unused(ptr);

TSGetDM(ts, );

DMCompositeGetEntries(pack, , NULL);

VecCopy(ydot, f);

DMCompositeGetLocalVectors(pack, , );
DMCompositeGetLocalVectors(pack, , );
DMCompositeScatter(pack, y, yda, yred);

VecSet(ydotda, ZERO);
VecSet(ydotred, ZERO);

func(t, yda, yred, ydotda, ydotred, da);

VecScale(ydotda, -ONE);
VecScale(ydotred, -ONE);

DMCompositeGather(pack, ADD_VALUES, f, ydotda, ydotred);
DMCompositeRestoreLocalVectors(pack, , );
DMCompositeRestoreLocalVectors(pack, , );

return 0;
}

void func(double t, Vec y, Vec aux, Vec ydot, Vec auxdot, DM da)
{
int rank;
PetscScalar ydotarr, *auxdotarr;
const PetscScalar yarr, *auxarr;
PetscInt i, j, k, is, js, ks, n, m, p, ys, ye;

unused(t);

MPI_Comm_rank(PETSC_COMM_WORLD, );

VecGetArrayRead(aux, );

DMDAGetCorners(da, , , , , , );

DMDAVecGetArrayDOFRead(da, y, );
DMDAVecGetArrayDOF(da, ydot, );

for (k = ks; k < ks + p; ++k)
{
for (j = js; j < js + m; ++j)
{
for (i = is; i < is + n; ++i)
{
double temp = auxarr[AUXTEMP];
double energy = 100.;
double R = NU * exp(-energy / temp);

ydotarr[k][j][i][0] += R;
ydotarr[k][j][i][1] -= R;
}
}
}

DMDAVecRestoreArrayDOFRead(da, y, );
DMDAVecRestoreArrayDOF(da, ydot, );

VecGetOwnershipRange(auxdot, , );
VecGetArray(auxdot, );

if (rank == 0)
{
for (i = ys; i < ye; ++i)
{
switch (i)
{
case AUXTEMP:
auxdotarr[i-ys] = HEATING;
break;
}
}
}

VecRestoreArray(auxdot, );
VecRestoreArrayRead(aux, );
}

int main(int argc, char *argv[])
{
DM da, red, pack;
Vec y, yda, yred;
PetscScalar yarr, *auxarr;
PetscInt i, j, k, is, js, ks, n, m, p, as, ae;
double t = ZERO, abstol = 1.e-6, reltol = 1.e-3;
int rank, tt = 0, ndof = 2, naux = 1, nx = 3, ny = 3, nz = 3;
PetscBool implicit = PETSC_TRUE;

TS ts;
TSAdapt adapt;
TSConvergedReason reason;

PetscInitialize(, , NULL, NULL);

Re: [petsc-users] PCMGSetGalerkin() new inputs

[Jed Brown via petsc-users]

This wasn't explained well in the commit message.  The old code used the
Galerkin procedure on the "Pmat" (preconditioning matrix; which may or
may not be the same as the Amat) and set the result as both Amat and
Pmat of the coarse grid.  The new code allows you to specify.  If your
Amat and Pmat are the same on the finest level, then BOTH will offer the
same behavior as before.

Myriam Peyrounette via petsc-users  writes:

> Hi,
>
> I am currently comparing two codes based on PETSc. The first one uses
> PETSC 3.6.4 and the other one PETSc 3.10.2.
>
> I am having a look at the use of the function PCMGSetGalerkin(). With
> PETSc 3.6, the input is a boolean, while it is either
> PC_MG_GALERKIN_MAT, PC_MG_GALERKIN_PMAT or PC_MG_GALERKIN_BOTH with PETSc 
> 3.10.
>
> Which one of these new inputs should I use to have the same
> configuration as with PETSc 3.6?
>
> Thanks in advance,
>
> -- 
> Myriam Peyrounette
> CNRS/IDRIS - HLST
> --

[petsc-users] PCMGSetGalerkin() new inputs

[Myriam Peyrounette via petsc-users]

Hi,

I am currently comparing two codes based on PETSc. The first one uses
PETSC 3.6.4 and the other one PETSc 3.10.2.

I am having a look at the use of the function PCMGSetGalerkin(). With
PETSc 3.6, the input is a boolean, while it is either
PC_MG_GALERKIN_MAT, PC_MG_GALERKIN_PMAT or PC_MG_GALERKIN_BOTH with PETSc 3.10.

Which one of these new inputs should I use to have the same
configuration as with PETSc 3.6?

Thanks in advance,

-- 
Myriam Peyrounette
CNRS/IDRIS - HLST
--




smime.p7s
Description: Signature cryptographique S/MIME

Re: [petsc-users] saving results

[Matthew Knepley via petsc-users]

On Wed, Feb 20, 2019 at 4:43 AM Sal Am  wrote:

> Hi Matthew you were right,
>
> The matrix I have is very ill conditioned and my supervisor gave it for
> testing purposes. Having said that, I was able to solve it previously
> however, for some reason it said convergence reached at e-3 even though the
> default convergence tol is e-5 which is why I put rel_tol to e-7 and it
> crashed.
>
> Just to get something solving so you can check the discretization, I would
>>
>>   a) Use SuperLU or MUMPS, -pc_type lu
>>
>>   b) Then to get bigger use ASM/LU, -pc_type asm -pc_sub_type lu
>>
>> Then finally
>>
>>   c) For very big problems, I suspect that PCPATCHcould be made to work
>> with the right choices, but this is a research problem
>>
>
> That is very interesting, I have tried a) before and it does not work
> (lack of memory presumably). I will try to do b), but I am very interested
> in c). the final BOSS fight so to say is the 30million finite element model
> of the JET A2 Antenna so it would be great if I can utilise PETSc solving
> it. (Perhaps I should open another ticket for that later on)
>

Yep. The reason I thought it might work is that there is theory for it
working on positive semi-definite matrices (matrices which are
almost nice but can have a large null space). The idea is to capture the
local nullspace in each patch. There is an excellent writeup
of the idea, with all appropriate references, here:
https://arxiv.org/abs/1810.03315

But to stay on topic I have only used VecView to actually view the solution
> *after* it was done solving. Is there an example on how people actually
> utilise VecView or Jed's suggestion of KSPMonitor to have checkpoints
> *during* iterations so we can pick up wherever it crashed?
>

Sure. Here we set a custom monitor


https://bitbucket.org/petsc/petsc/src/6f3214441ec7b7dca9916d2bd112d01eef3e7185/src/ksp/ksp/examples/tutorials/ex9.c#lines-111

Just stick your VecView in there. In fact, I would use VecViewFromOptions()
so you can decide what viewer to use on the command line.

  Thanks,

 Matt


> Thanks
>
> On Tue, Feb 19, 2019 at 1:49 PM Matthew Knepley  wrote:
>
>> On Tue, Feb 19, 2019 at 8:21 AM Sal Am  wrote:
>>
>>> Matthew maybe indeed I was not clear and wrong in the wording. What I
>>> meant is the matrix is undefined, but since it is a finite element method
>>> (using second order elements). The pattern is like this:
>>>
>>
>> Okay, the paper says that the formulation of Maxwell in frequency space
>> is "well conditioned", and in the paper they solve using
>> BiCG/Jacobi:
>>
>>   -ksp_type bicg -pc_type jacobi
>>
>> Since you are taking thousands of iterates:
>>
>>   a) They are mistaken and the formulation is not well-conditioned
>>
>>   b) You have a bug in the discretization
>>
>>   c) You are outside the parameter range/boundary conditions/forcing they
>> were talking about for conditioning
>>
>> Just to get something solving so you can check the discretization, I would
>>
>>   a) Use SuperLU or MUMPS, -pc_type lu
>>
>>   b) Then to get bigger use ASM/LU, -pc_type asm -pc_sub_type lu
>>
>> Then finally
>>
>>   c) For very big problems, I suspect that PCPATCHcould be made to work
>> with the right choices, but this is a research problem
>>
>>   Thanks,
>>
>>  Matt
>>
>>
>>> [image: image.png]
>>> We are indeed solving Maxwell's equations. I have added the paper.
>>>
>>> On Tue, Feb 19, 2019 at 12:15 PM Matthew Knepley 
>>> wrote:
>>>
 On Tue, Feb 19, 2019 at 4:12 AM Sal Am via petsc-users <
 petsc-users@mcs.anl.gov> wrote:

> It is a finite-element problem of an RF antenna dielectric interaction
> where all the non-zero elements are on the diagonal of the sparse matrix
> (if that is relevant).
>

 I am unsure whether this is what you mean. If you matrix is diagonal,
 you can invert it exactly in a single step so you
 do not need a Krylov solver. What equations are you solving? Maxwell?
 Electrostatic? What finite element are you using?
   Thanks,

 Matt


> More about the matrix: 25Mx25M, 2B non-zeros. I checked the
> KSPMonitor, it seems that I have to write my own routine?
>
>> Running for a Krylov method for tens of thousands of iterations is
>> very rarely recommended.
>>
> GAMG and BCGS is the only ones that have actually worked for me so
> far, I increased the max_it because it was not enough with the default 
> one.
> With the default tolerance I got ||r(i)||/||b|| in the order of 10^-3, but
> I need more accuracy so I increased the tol to 10^-7 and that is when it
> crashed after ~51000 iterations.
>
> @Barry
>
>> I'm guessing you mange your own time stepping (and nonlinear solver
>> if there is one)?
>>
>
> It is the default from the solver (attached the short code).
>
>> As Jed says it doesn't make sense to save "partial" solutions within
>> the linear solver.  Just save it