[Nek5000-users] Interfacing with Arpack

[nek5000-users at lists.mcs.anl.gov]

Hi Giuseppe,

You're close.

The key point to understand is that application of Q^T arises from continuity
requirements on the test functions (that we normally call "v" in an abstract setting).
Everything in Nek is cast in the weak form:   Find u \in X^N_0 such that, for all v \in X^N_0,

        a( v, u )_N =   a( v, u* )_N

where u* is the exact solution and a( . , . )_N is the inner-product of interest, with discrete
quadrature.   For example, if 

     a( v , u ) := \int_Omega  grad v . grad u  dV          (1)

then a( v , u )_N would be the same bilinear form with integration replaced by numerical
quadrature.    Note that an L2 projection would be of the form:


      ( v , u ) = ( v , u* )

where I've dropped the _N on the inner product on the working assumption that we
always use quadrature.

Since u,v are continuous in Nek, the bilinear forms read:

           v^T  Q^T M_L Q u =   v^T Q^T M_L  u*_L

and, in nek,

             u = Binvm1 * QQ^T  BM1  u*

Note that (1), above, has the mass matrix build in because of the volume
integral.    So, technically, if you are doing dssum, you should have a mass
matrix present, because dssum relates to projection against the test functions
v, which implies a volume-weighted average of the equations.

Sorry if this isn't very clear...

Paul




________________________________________
From: nek5000-users-bounces at lists.mcs.anl.gov [nek5000-users-bounces at lists.mcs.anl.gov] on behalf of nek5000-users at lists.mcs.anl.gov [nek5000-users at lists.mcs.anl.gov]
Sent: Monday, October 20, 2014 8:21 AM
To: nek5000-users at lists.mcs.anl.gov
Subject: Re: [Nek5000-users] Interfacing with Arpack

Dear Paul,
thank you for the answer.
If my understanding is correct, since each matrix-vector product (coded
as matrix-matrix product) is done on the local variables, I should call
dssum after each matrix-vector operation.
For instance, if I want to apply the Helmholtz matrix to a vector u and
store the result in v, I should do:
call axhelm(v,u,h1,h2,imsh,1)
call dssum(v,nx1,ny1,nz1)
is it correct?
If I want to apply the gradient to a vector (always for eigenvalue
computation purposes) shoud I do something like:
       call gradm1(dudx,dudy,dudz,u)
       call dssum(dudx,nx1,ny1,nz1)
       call dssum(dudy,nx1,ny1,nz1)
       call dssum(dudz,nx1,ny1,nz1)
       bufx = glsc2(u0,dudx,n)           !bufx = u0.du/dx
or am I missing something?

Giuseppe



On 20/10/2014 12:58, nek5000-users at lists.mcs.anl.gov wrote:
> Giuseppe,
>
> The key thing to understand about Nek is that variables are _always_ stored in their local format,
> which implies you have redundant values at the interface such that they are ready for parallel operator
> evaluation.
>
> Let u_L denote such a vector with redundant values (L for Local) and let u be the corresponding
> vector without the redundantly stored data.  Nominally, there is a matrix Q that maps u to u_L:
>
> (1)       u_L = Q u
>
> and the matrix-vector product involving the mass matrix, of the form w = M u, can be expressed as:
>
> (2)       w = M u =  Q^T M_L Q u  = Q^T  M_L u_L
>
> where M_L = block_diag(M^e) , e=1,...,Nelv is the block-diagonal matrix consisting of local mass
> matrices within each element, e.   In Nek, M^e is also block diagonal, so the whole _unassembled_
> stiffness matrix, M_L, is diagonal.   It is this matrix that is stored in bm1  ("B" on Mesh 1, the velocity/temp
> mesh).
>
> Notice that, as written, the output of (2) violates our "store only local vectors" policy, so we invoke
> (1) to correct this:
>
>
> (3)    w_L = Q w = Q M u =  Q Q^T  M_L u_L.
>
>
> We refer to QQ^T as direct-stiffness summation  (dssum, in nek).   It is the operation that
> glues things together.
>
> So, suppose you wish to find a continuous function that is close to a discontinuous function
> in nek.   Let f be the continuous thing and g be the discontinuous one.   You can do this via:
>
> (4)      f_L  = Q  M^-1   Q M_L g_L.
>
> In Nek, this would be written as:
>
>            n = nx1*ny1*nz1*nelv
>
>            call col3(f,g,bm1,n)
>
>            call dssum(f,nx1,ny1,nz1)
>
>            call col2    (f,binvm1,n)
>
> Here, binvm1 is  effectively   binvm1 =  1./ (QQT*bm1)   , if we view bm1 and binvm1 as
> vectors, rather than diagonal matrices.
>
> Paul
>
>
>
>
>
> ________________________________________
> From: nek5000-users-bounces at lists.mcs.anl.gov [nek5000-users-bounces at lists.mcs.anl.gov] on behalf of nek5000-users at lists.mcs.anl.gov [nek5000-users at lists.mcs.anl.gov]
> Sent: Monday, October 20, 2014 2:03 AM
> To: nek5000-users at lists.mcs.anl.gov
> Subject: [Nek5000-users] Interfacing with Arpack
>
> Dear users and developers,
> I would like to solve a generalized eigenvalue problem with Arpack,
> based on Nek subroutines. Using reverse communication interface, I only
> need to express the action of the matrices on the vectors, but I am not
> sure on how to apply the mass matrix and the inverse mass matrix in the
> Pn-Pn formulation.
> Is it sufficient to write:
> call col2(p1,bm1,n)      !apply mass matrix to vector p1
> ...
> call invcol2(p1,bm1,n)  ! apply inverse mass matrix to vector p1
> or should I do:
> call rzero(h1,n)
> call rone(h2,n)
> call hmholtz('vel',p,p1,h1,h2,pmask,vmult,imsh,tolh,maxiter,isd)
> to apply inverse mass matrix to p1, or should I do something else?
> Also, from a related question from a couple of years ago, should I call
> col2(v,vmask,n) after each matrix-vector product?
> Thank you in advance.
> Best,
> Giuseppe
> _______________________________________________
> Nek5000-users mailing list
> Nek5000-users at lists.mcs.anl.gov
> https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
> _______________________________________________
> Nek5000-users mailing list
> Nek5000-users at lists.mcs.anl.gov
> https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users

_______________________________________________
Nek5000-users mailing list
Nek5000-users at lists.mcs.anl.gov
https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users