[petsc-users] Mass matrix with PetscFE

[Julian Andrej]

Hey,

yes, this issue was resolved a few weeks after the mail. I just tried
again after some DMPlex commits ;).

Thanks!

On Thu, Jul 21, 2016 at 3:04 PM, Matthew Knepley <knepley at gmail.com> wrote:
> On Mon, Mar 7, 2016 at 6:21 PM, Julian Andrej <juan at tf.uni-kiel.de> wrote:
>>
>> Any news about this? I've seen you merged the dmforest branch into next.
>
>
> I am going through my mail, and see that this might have been dropped. Has
> your
> problem been solved?
>
>   Sorry for the delay,
>
>      Matt
>
>>
>>
>> On 2016-02-26 01:22, Matthew Knepley wrote:
>>>
>>> I am sorry about the delay. I have your example working but it exposed
>>> a bug in Plex so I need to push the fix first. I should have
>>> everything for you early next week.
>>>
>>>   Thanks
>>>
>>>      Matt
>>>
>>> On Feb 25, 2016 2:04 AM, "Julian Andrej" <juan at tf.uni-kiel.de> wrote:
>>>
>>>> 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 [1]: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 [2]: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 [3] <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
>>>
>>>
>>>
>>> Links:
>>> ------
>>> [1] tel:24.02.2016%2020
>>> [2] tel:24.02.2016%2014
>>> [3] tel:24.02.2016%2013
>
>
>
>
> --
> 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