[Nek5000-users] Element connectivity information

Thu Mar 15 10:53:05 CDT 2012

Dear Outi,

My suspicion is that there is an issue with the way Arnoldi
is viewing the matrix-vector product.  A velocity vector that
is discontinuous is inadmissable.  In the example below, that
would correspond to failing to multiply by P, which would mean
that you would be working with a matrix that is different from
A'.

A proper matvec in Nek is of the form:

       call axhelm(v,u,...)
       call dssum (v...)
       call col2  (v,mask,...)

The dssum() ensures that v is in H1.  The mask ensures that it is
in H^1_0.   Note that u should also be in H^1_0 if it is to be
a legitimate element of the Krylov subspace in which the Arnoldi
scheme is operating.

(Obviously, a matvec for your case is more complex... but I just
wanted to show the essential ingredients for the simple case of
the Helmholtz operator.)

Best regards,

Paul

On Thu, 15 Mar 2012, nek5000-users at lists.mcs.anl.gov wrote:

> Dear Paul,
>
> Thank you for the input. I will think about your answers and return to you.
>
> I should mention that it is true that the zero eigenvalues have eigenvectors 
> that are discontinuous at the element boundaries.
> My idea was that velocity vectors that are slightly discontinuous across the 
> element boundaries might be returned by Arnoldi, if not else, due to finite 
> precision arithmetics. (And since Arnoldi is looking for the unstable or 
> least stable eigenvalues, zero eigenvalues would often be returned among the 
> first ones.)
>
> But I understand that your experience tells you the problem is somewhere 
> else.
>
> Outi
>
> 15 Mar 2012 kl. 13:44 skrev nek5000-users at lists.mcs.anl.gov:
>
>> 
>> Dear Outi,
>> 
>> Just a follow-up on my earlier comments about Arnoldi.
>> 
>> With iterative methods in general, there is absolutely no
>> need to condense the redundant nodes because the solver
>> has no concept of the dimension of the system ("n" is not
>> an explicit part of the algorithm).   The solver only knows about the rank 
>> that it discovers during the orthogonalization
>> process.  So, as long as every input vector in the SEM case is continuous, 
>> arpack will not distinguish whether you have the
>> redundant storage or not---it has no means for making this
>> discovery.
>> 
>> Hence, if you are seeing spurious modes, they are a result
>> of breaking the condition that the vectors are continuous
>> or of some other mode that is not being controlled.
>> A common source for such breakage is failure to dssum in
>> the appropriate place, particularly in a preconditioning
>> step, or somehow failing to propagate needed information.
>> My suspicion is that there is something else going on
>> that we should pay attention to.
>> 
>> Just to make this concrete, arpack will discover the eigenvalues of the 
>> following to matrices (of rank 2)
>> in exactly the same iteration count:
>> 
>>
>>    A(2x2)
>> 
>> and
>>
>>    A' := PAP^T
>> 
>> 
>> where
>>
>>        / 1 0 \
>>    P = | 0 1 |
>>        \ 0 1 /
>> 
>> note that the latter case makes vectors v' = P v have
>> redundant 2nd and 3rd entries.
>> 
>> This being said, you have to pay attention to any preconditioning steps 
>> (particularly any involving direct methods, where knowing "n" is 
>> important), to
>> ensure that each vector in the iteration step is in the
>> column space of P.
>> 
>> (In nek, the equivalent of "P" is the matrix that does
>> the global to local map.... and in fact, there is never a matrix "A" of the 
>> type described above.)
>> 
>> Note that it might appear that A' has a spurious
>> eigenmode, which is true, but it is in fact not
>> allowed because that eigenmode is not in the range
>> of P (i.e., would correspond to a discontinuous
>> velocity field, which would never be generated in
>> a proper setup).   Also note that there might be an
>> issue if you told Arpack that the rank of your matrix
>> is 3 instead of 2 (e.g., it would look for 3 eigenvalues).
>> However, since you are only looking for m << n eigenvalues,
>> there is no reason to suspect that this is an issue.
>> 
>> I hope this helps.
>> 
>> Cheers,
>> 
>> Paul
>> 
>> On Thu, 15 Mar 2012, nek5000-users at lists.mcs.anl.gov wrote:
>> 
>>> Dear Paul,
>>> 
>>> I am using the parallel version of ARPACK. It is currently only called 
>>> inside userchk. When the Arnoldi operations finish, the Nek5000 velocity 
>>> vector is updated, and I force Nek5000 to start with the new vector as an 
>>> initial condition.
>>> (I am aware of that this is not an elegant implementation, but we thought 
>>> other options would imply modifying the source).
>>> 
>>> I have verified the code against the local Orr-Sommerfeld spectrum in a 
>>> channel flow (using periodic BC in the streamwise direction), and the wake 
>>> behind a cylinder.
>>> It seems like the spectrum is reproduced correctly, except from a few 
>>> spurious eigenvalues, that are always located "exactly" at zero (real and 
>>> imaginary part). Those do not appear for all cases I have studied, but for 
>>> many cases. So far, I have ignored the spurious eigenvalues, though they 
>>> make the convergence slower.
>>> 
>>> Best regards,
>>> Outi
>>> 
>>> 14 Mar 2012 kl. 22:23 skrev nek5000-users at lists.mcs.anl.gov:
>>> 
>>>> Dear Outi,
>>>> Can you tell me a bit more about the Arnoldi solver? (I know
>>>> about the method.  Here, I mean the specific source code that
>>>> you are using, what the interface to the code is, etc.)
>>>> Best regards,
>>>> Paul
>>>> On Wed, 14 Mar 2012, nek5000-users at lists.mcs.anl.gov wrote:
>>>>> Hi Paul,
>>>>> Thank you for a quick response.
>>>>> My current understanding is that the error occurs since Arnoldi treats 
>>>>> the velocities in the duplicated points as separate degrees of freedom 
>>>>> (which they are not in Nek5000, since they come out always equal).
>>>>> So, therefore I am not so sure that a weighting operation would remove 
>>>>> the problem.
>>>>> Anyway, I will test what you suggest in a few cases and tell you how it 
>>>>> comes out.
>>>>> But if it turns out that counting the points once is the only option 
>>>>> (asymmetrically, as you point out), is it possible to do this?
>>>>> (Another option is always that I make changes inside the Arnoldi 
>>>>> algorithm, but this might be a bit more elaborate.)
>>>>> Thanks a lot,
>>>>> Outi
>>>>> 14 Mar 2012 kl. 13:40 skrev nek5000-users at lists.mcs.anl.gov:
>>>>>> Hi Outi,
>>>>>> We encounter this problem all the time in conjugate gradient
>>>>>> and gmres iteration.
>>>>>> For reasons of symmetry (and, hence, commutativity and associativity,
>>>>>> which are two essential ingredients for parallel execution), I would
>>>>>> suggest that you not try to store the data only once, but rather
>>>>>> scale by the inverse counting-matrix (aka, multiplicity), so that 
>>>>>> interface points are not counted twice when computing your inner
>>>>>> products.
>>>>>> The inverse counting-matrix is stored in vmult() and tmult().  A
>>>>>> typical application is:
>>>>>> 
>>>>>> alpha = glsc3(vx,vmult,vx,n)   (n=nx1*ny1*nz1*nelv)
>>>>>> which gives the l2 inner product:  alpha := sum_i (vx_i)^2 for a vector 
>>>>>> VX.
>>>>>> Note that if you instead want the L2 inner product
>>>>>> 
>>>>>> alpha = \int_Omega (vx)^2 dV
>>>>>> you would compute this as:
>>>>>> 
>>>>>> alpha = glsc3(vx,bm1,vx,n)
>>>>>>
>>>>>>       = vx^T B vx
>>>>>> which is the inner product of vx with the (diagonal) mass matrix B.
>>>>>> If you want the average of the physical rms (say), it would be:
>>>>>> 
>>>>>> alpha = glsc3(vx,bm1,vx,n)
>>>>>> if (alpha.gt.0) alpha = sqrt(alpha/volvm1)
>>>>>> etc.
>>>>>> In summary, for linear algebra, use vmult(). For physics, use bm1()
>>>>>> (B, on Mesh 1).
>>>>>> Please let me know if this help to resolve your problem.
>>>>>> Best regards,
>>>>>> Paul
>>>>>> On Wed, 14 Mar 2012, nek5000-users at lists.mcs.anl.gov wrote:
>>>>>>> Hi Neks,
>>>>>>> I am using Nek5000 coupled to an Arnoldi algorithm for hydrodynamic 
>>>>>>> stability computations. In these, I form a vector of the velocities, 
>>>>>>> which I regularly pass to the Arnoldi solver. Then, a new velocity 
>>>>>>> vector is formed which is passed back to Nek5000.
>>>>>>> Now, the code works as it is, but in addition to the physical modes it 
>>>>>>> produces some spurious modes, which I think originate from the fact 
>>>>>>> that the velocities at the element-common boundaries are accounted for 
>>>>>>> twice. To fix this, I need to access the element connectivity 
>>>>>>> information. Can you please guide me to in which variable this is 
>>>>>>> stored?
>>>>>>> In particular, this is what I would like to do:
>>>>>>> -  add the element-common points to the local vector only for one of 
>>>>>>> two adjacent elements (and for a corner point, only one of the four 
>>>>>>> elements adjacent to it).
>>>>>>> - for the vector that Arnoldi returns, make sure that the new values 
>>>>>>> of the velocity are inserted also for the other adjacent element(s).
>>>>>>> Thanks a lot in advance!
>>>>>>> Outi
>>>>>>> _______________________________________________
>>>>>>> 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