[petsc-users] Mass matrix with PetscFE

[Julian Andrej]

After a bit of rethinking the problem, the discrepancy between the size 
of matrix A and the mass matrix M arises because of the Dirichlet 
boundary conditions. So why aren't the BCs not imposed on the mass 
matrix? Do I need to handle Dirichlet BCs differently in this context 
(like zero rows and put one the diagonal?)

On 24.02.2016 20:54, juan wrote:
> I attached another example which creates the correct mass matrix
> but also overwrites the DM for the SNES solve. Somehow i cannot manage
> to really copy the DM to dm_mass and use that. If i try to do that with
> DMClone(dm, &dm_mass) i get a smaller mass matrix (which is not of size A).
>
> Maybe this helps in the discussion.
>
> Relevant code starts at line 455.
>
> On 2016-02-24 15:03, Julian Andrej wrote:
>> Thanks Matt,
>>
>> I attached the modified example.
>>
>> the corresponding code (and only changes to ex12) is starting at line
>> 832.
>>
>> It also seems that the mass matrix is of size 169x169 and the
>> stiffness matrix is of dimension 225x225. I'd assume that if i
>> multiply test and trial function i'd get a matrix of same size (if the
>> space/quadrature is the same for the stiffness matrix)
>>
>> On 24.02.2016 14:56, Matthew Knepley wrote:
>>> On Wed, Feb 24, 2016 at 7:47 AM, Julian Andrej <juan at tf.uni-kiel.de
>>> <mailto:juan at tf.uni-kiel.de>> wrote:
>>>
>>>     I'm now using the petsc git master branch.
>>>
>>>     I tried adding my code to the ex12
>>>
>>>        DM dm_mass;
>>>        PetscDS prob_mass;
>>>        PetscFE fe;
>>>        Mat M;
>>>        PetscFECreateDefault(dm, user.dim, 1, PETSC_TRUE, NULL, -1, &fe);
>>>
>>>        DMClone(dm, &dm_mass);
>>>        DMGetDS(dm_mass, &prob_mass);
>>>        PetscDSSetDiscretization(prob_mass, 0, (PetscObject) fe);
>>>        PetscDSSetJacobian(prob_mass, 0, 0, mass_kernel, NULL, NULL,
>>> NULL);
>>>        DMCreateMatrix(dm_mass, &M);
>>>
>>>        MatSetOptionsPrefix(M, "M_";)
>>>
>>>     and receive the error on running
>>>     ./exe -interpolate -refinement_limit 0.0125 -petscspace_order 2
>>>     -M_mat_view binary
>>>
>>>     WARNING! There are options you set that were not used!
>>>     WARNING! could be spelling mistake, etc!
>>>     Option left: name:-M_mat_view value: binary
>>>
>>>     I don't know if the matrix is actually there and assembled or if the
>>>     option is ommitted because something is wrong.
>>>
>>>
>>> Its difficult to know when I cannot see the whole code. You can always
>>> insert
>>>
>>>    MatViewFromOptions(M, NULL, "-mat_view");
>>>
>>>     Using
>>>     MatView(M, PETSC_VIEWER_STDOUT_WORLD);
>>>
>>>     gives me a reasonable output to stdout.
>>>
>>>
>>> Good.
>>>
>>>     But saving the matrix and analysing it in matlab, results in an all
>>>     zero matrix.
>>>
>>>     PetscViewerBinaryOpen(PETSC_COMM_WORLD, "Mout",FILE_MODE_WRITE,
>>>     &viewer);
>>>     MatView(M, viewer);
>>>
>>>
>>> I cannot explain this, but it has to be something like you are viewing
>>> the matrix before it is
>>> actually assembled. Feel free to send the code. It sounds like it is
>>> mostly working.
>>>
>>>    Matt
>>>
>>>     Any hints?
>>>
>>>
>>>     On 24.02.2016 13 <tel:24.02.2016%2013>:58, Matthew Knepley wrote:
>>>
>>>         On Wed, Feb 24, 2016 at 6:47 AM, Julian Andrej
>>>         <juan at tf.uni-kiel.de <mailto:juan at tf.uni-kiel.de>
>>>         <mailto:juan at tf.uni-kiel.de <mailto:juan at tf.uni-kiel.de>>>
>>> wrote:
>>>
>>>              Hi,
>>>
>>>              i'm trying to assemble a mass matrix with the
>>> PetscFE/DMPlex
>>>              interface. I found something in the examples of TAO
>>>
>>> https://bitbucket.org/petsc/petsc/src/da8116b0e8d067e39fd79740a8a864b0fe207998/src/tao/examples/tutorials/ex3.c?at=master&fileviewer=file-view-default
>>>
>>>
>>>              but using the lines
>>>
>>>              DMClone(dm, &dm_mass);
>>>              DMSetNumFields(dm_mass, 1);
>>>              DMPlexCopyCoordinates(dm, dm_mass);
>>>              DMGetDS(dm_mass, &prob_mass);
>>>              PetscDSSetJacobian(prob_mass, 0, 0, mass_kernel, NULL,
>>>         NULL, NULL);
>>>              PetscDSSetDiscretization(prob_mass, 0, (PetscObject) fe);
>>>              DMPlexSNESComputeJacobianFEM(dm_mass, u, M, M, NULL);
>>>              DMCreateMatrix(dm_mass, &M);
>>>
>>>              leads to errors in DMPlexSNESComputeJacobianFEM (u is a
>>>         global vector).
>>>
>>>              I don't can understand the necessary commands until
>>>              DMPlexSNESComputeJacobianFEM. What does it do and why is it
>>>              necessary? (especially why does the naming involve SNES?)
>>>
>>>              Is there another/easier/better way to create a mass
>>> matrix (the
>>>              inner product of the function space and the test space)?
>>>
>>>
>>>         1) That example needs updating. First, look at SNES ex12 which
>>>         is up to
>>>         date.
>>>
>>>         2) I assume you are using 3.6. If you use the development
>>>         version, you
>>>         can remove DMPlexCopyCoordinates().
>>>
>>>         3) You need to create the matrix BEFORE calling the assembly
>>>
>>>         4) Always always always send the entire error messge
>>>
>>>             Matt
>>>
>>>              Regards
>>>              Julian Andrej
>>>
>>>
>>>
>>>
>>>         --
>>>         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
>>>
>>>
>>>
>>>
>>> --
>>> 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