[Nek5000-users] Interfacing with Arpack

nek5000-users at lists.mcs.anl.gov nek5000-users at lists.mcs.anl.gov
Mon Oct 27 11:17:47 CDT 2014 

Hi Adam,
thanks for sharing the code, I will study it, and if I can produce 
something useful I will contact you.
Best regards,
Giuseppe



On 27/10/2014 16:47, nek5000-users at lists.mcs.anl.gov wrote:
> Hi Guiseppe
>
> I've just noticed you work on eigenvalue problem with arpack for 
> nek5000 p_n-p_n. We've done similar work for p_n-p_n-2 and we are very 
> interested in your work and we would like to share the code if it is 
> possible. Our implementation of ext_cyl example you can find at
> ftp://ftp.mech.kth.se/pub/adam/Nek5000/toolbox/
> Best regards
> Adam
>
> On 22/10/14 13:31, nek5000-users at lists.mcs.anl.gov wrote:
>>
>> 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
>>>
>> _______________________________________________
>> 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

More information about the Nek5000-users mailing list