[Nek5000-users] Interfacing with Arpack
nek5000-users at lists.mcs.anl.gov
nek5000-users at lists.mcs.anl.gov
Mon Oct 27 10:47:10 CDT 2014
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
More information about the Nek5000-users
mailing list