[Nek5000-users] Interfacing with Arpack

Wed Oct 22 01:53:57 CDT 2014

Hi Paul, thanks again for your reply. I think I am close to a solution 
now. One more thing: I think I need to mask out Dirichlet bc nodes when 
applying the matrix-vector products, right? How can I do this? Is it 
sufficient to call col2(u,vmask,n) after the matrix-vector product?
Best,
Giuseppe

On 20/10/2014 16:37, nek5000-users at lists.mcs.anl.gov wrote:
> 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
> _______________________________________________
> Nek5000-users mailing list
> Nek5000-users at lists.mcs.anl.gov
> https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users