[Nek5000-users] Element connectivity information

Thu Mar 15 08:44:20 CDT 2012

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