[Nek5000-users] Interfacing with Arpack

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

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