[Nek5000-users] Element connectivity information

Thu Mar 15 09:16:00 CDT 2012

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