[Nek5000-users] Interfacing with Arpack

Wed Oct 22 06:31:08 CDT 2014

Giuseppe,

There are several masks:  tmask()

v1mask for vx
v2mask for vy
v3mask for vz

pmask for pressure

tmask has several fields in the situation where
you have different boundary conditions on different
passive scalars.

Note that repeated application of one of the SEM operators
requires injection into H^1, which is effected through dssum
+ mass matrix, as discussed earlier.   Subroutine ALPHAM1
in eigsolv illustrates this.  (It's an older routine, so it
uses the multiplicity array when computing this inner products.
That would not be necessary, however, if the inner products
had been computed prior to dssum.)

Paul

On Wed, 22 Oct 2014, nek5000-users at lists.mcs.anl.gov wrote:

> 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
>
> _______________________________________________
> Nek5000-users mailing list
> Nek5000-users at lists.mcs.anl.gov
> https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
>