On Tue, Sep 1, 2015 at 12:35 PM, Julian Andrej <j...@tf.uni-kiel.de> wrote:
> On Tue, Sep 1, 2015 at 6:15 PM, John Peterson <jwpeter...@gmail.com>
> wrote:
>
> >
> >
> > On Tue, Sep 1, 2015 at 8:29 AM, Julian Andrej <j...@tf.uni-kiel.de>
> wrote:
> >
> >> I'm trying to get into libMesh, so i tried to convert an existing
> project.
> >>
> >> So far i could reproduce the lid driven cavity case for the Stokes
> >> equations like in
> >> example examples/systems_of_equations/systems_of_equations_ex1.
> >>
> >> Now i'm trying to build a mass matrix so i can solve the generalized
> >> eigenproblem in order to retrieve the physical eigenvalues in the
> >> following
> >> domain:
> >>
> >> build_square (mesh,
> >> 16, 16,
> >>
> >> -1., 1.,
> >>
> >> -1., 1.,
> >>
> >> TRI6);
> >>
> >> So the eigenvalues should be (rounded to int) [13, 23, 23, 32, 38,
> 41...].
> >>
> >> My approach is that i add a matrix
> >>
> >> system.add_matrix("mass");
> >>
> >> in the assemble loop i add:
> >>
> >> DenseMatrix<Number> Me;
> >> Me.resize (n_dofs, n_dofs);
> >>
> >> and to the integration parts:
> >> Me(i,j) += JxW[qp]*(phi[i][qp]*phi[j][qp]);
> >>
> >> then i close the matrices
> >> system.get_matrix("mass").close();
> >>
> >> and convert them to PETSc so i can use my own SLEPc object
> >> Mat A = (static_cast<PetscMatrix<Number> *>(system.matrix))->mat();
> >>
> >> Mat M = (static_cast<PetscMatrix<Number>
> >> &>(system.get_matrix("mass"))).mat();
> >>
> >> The problem is, the calculated eigenvalues are not the ones i expect (so
> >> they're wrong).
> >> This method worked for me for the poisson equation.
> >> Do i miss something?
> >>
> >
> > So are you trying to solve the generalized eigenvalue problem
> >
> > A*x = lambda*B*x,
> >
>
> > where A corresponds to the Stokes operator and B to the "mass" matrix,
> > which should have a zero block for the pressure equation?
> >
>
> Exactly
>
>
> > If yes, then
> > .) Are you building both matrices, since you only talk about the mass
> > matrix?
> >
>
> Yes i build the Stiffness matrix based on the example code
>
>
> > .) Are you building the correct mass matrix for Stokes?
> >
>
> That is the main question i guess. I'm adding the code
>
> Me(i,j) += JxW[qp]*(phi[i][qp]*phi[j][qp]);
>
> to the Kuu and Kvv loops. I've also tried adding it to the Kup, Kvp, Kpv,
> Kvp, as well but the quadrature seems to be zero there (using mat view from
> petsc and viewing the matrix structure).
>
>
> >
> > Another question is BCs: what are they, and how are you imposing them
> > (e.g. penalty method, constraints, etc.)? I can't remember if BCs should
> > be imposed in both the A and B matrices or not...
> >
>
> I'm imposing pure zero dirichlet using the penalty method as in the example
> code. I tried applying them just to A and also tried A and M.
>
> You can take a look at what i've got in that paste
> http://hastebin.com/raw/ekasibazan. I don't know if it is desired that i'd
> attach the code (i could do so for further reference). It is basically the
> example code with a few changes.
>
Regarding BCs: Constraint rows introduce spurious eigenvalues. For
Dirichlet boundary conditions, I would recommend following the approach in
eigenproblems_ex3 which uses a CondensedEigenSystem, which eliminates
Dirichlet dofs from the system before doing the solve.
David
------------------------------------------------------------------------------
_______________________________________________
Libmesh-users mailing list
Libmesh-users@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/libmesh-users
