[Nek5000-users] Problem with intp_rstd

[nek5000-users at lists.mcs.anl.gov]

Dear Juan,

The usage model for intp_rstd() is to map back and forth between a std grid and a fine grid,
in the context of Galerkin projection.  Here is the basic idea:

If you want to find u \in PN such that    ||(u*-u)|| \leq ||(u*-v)|| for all v \in PN, that is, such that u
is the closest element in PN to u*, then you have to choose a norm.

Let's suppose you are choosing the L2 norm.   Let's also assume u*=vx^2, which is a polynomial of
degree 2N.

If u is the closest element of PN to u*, it means that  (u*-u) \perp PN, which is to say,

      (v, (u-u*) ) = 0 for all v in PN   [ i.e.,  (v,e)=0, where e=u-u* is the error and (.,.) is the L2 inner-product].

Assuming there are n basis functions in PN, we have:

          u = sum_j=1^n  u_j  phi_j (X)

and thus we have n unknown basis coefficients, u_j, j=1,...n.

Moreover, the projection (i.e., best-fit) statement gives us n equations:

       (phi_i , u-u* ) = 0 ,   i=1,...,n

or,

       (phi_i, u )   = (phi_i , u*)   , i=1,...,n

Bescause the inner product  (f,g) := \int_Omega f g dX is linear, we  can right the term on the left as

     sum_j  (phi_i, phi_j )  u  =  (phi_i , u*)  =:  b_i ,   i=1,...,n

so   B u = b

where B_{ij} :=  (phi_,phi_j) := \int_Omega phi_i , phi_j dX.

Thus, we have two types of integrals to compute, B_{ij} and b_i.   Let's take care of b_i first:

   b_i = \int_Omega v u* dX

If u* \in P_2N and v \in P_N, then we need a quadrature rule that is exact for all polynomials of
degree M \leq 3N.    Gauss (not Gauss Lobatto) quadrature on (3M/2+1) points is exact for all
polynomials of degree  2(3M/2+1) - 1 = 3M+1, so that will suffice.   To compute b_i, we thus
do:

                b_i = sum_k=1^m  phi_i (\xi_k)  rho_k  vx(\xi_k) vx(\xi_k)

where \xi_k are the Gauss quadrature points in the domain, m=number of Gauss points,
rho_k are the Mth-order Gauss quadrature weights.

Note that, it is important to interpolate vx to the \xi_k, and not u*, because information is
otherwise already lost.    Note that if

      vx(X) = sum_j  \phi_j(X) vx_j

then

      vx(\xi_k) = \sum \phi_j(\xi_k) vx_j  =     J vx

where J_{kj}  :=  \phi_j(\xi_k) is the interpolation matrix from the Nth-order Gauss Lobatto (GL)
points to the Mth-order Gauss (G) quadrature points.    Careful inspection of of the above expression
shows that the rhs ( b = {b_i} ) is given by:

          b = J^T B_M V  (J vx)

where V = diag ( vx(\xi_k) ) is just a diagonal matrix with vx on the fine quadrature mesh.

So, ignoring momentarily the diagonal matrix V, the process of generating the right hand side,
with entries on the Nth-order Gauss points is, in words:

          Interpolate vx onto fine mesh:   vx -->  J vx

          Collocate with V and the diagonal matrix of quadrature weights, B_M

          Map back to the Nth-order GL mesh:       J^T ( B_M V (J vx ) )

OK.   So,

      intp_rstd(...,0) applies J.

      intp_rstd(...,1) applies J^T

Now, you still haven't finished the projection process.

To get u, you must do the following:

      call col3(u,binvm1,rhs,n)

which applies B^{-1} to get u  (from, as it were, u*= vx * vx ).


As a final finishing touch, the "col3" statement will be exact if the product u*v is in P_2N-1.  Unfortunately,
it is in P_2N, because u and v are both polynomials of degree N.

You have the option of computing the non diagonal matrix and solving that system via preconditioned
conjugate gradients, or of just allowing that small discrepancy.   If you've come this far, I would
recommend the former.   I have code that can do that for you.

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: Sunday, April 09, 2017 2:12 AM
To: nek5000-users at lists.mcs.anl.gov
Subject: [Nek5000-users] Problem with intp_rstd


Dear all,


I've been playing around with the 'intp_rstd' subroutine to do calculations on the dealiasing grid (lxd*lxy*lzd). I can map the velocity field onto the fine grid (vx --> vxd) and when I print out values it seems fine. Then I try to calculate vxd**2, and still seems fine. However, when I try to map back to the original grid, the results are complete garbage. Here is a snippet of my code.


        do e=1,nelv    ! mapping to fine grid
          call intp_rstd(vxd(1,1,1,e),vx(1,1,1,e),nx1,nxd,if3d,0) ! 0 --> forward
        enddo

          call vsq(vxd,ntotd)  ! calculating vxd**2

        do e=1,nelv    ! mapping back to original grid
          call intp_rstd(vx(1,1,1,e),vxd(1,1,1,e),nx1,nxd,if3d,1) ! 0 --> backwards
        enddo


Am I doing something wrong?


Thanks,


Juan Diego
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/nek5000-users/attachments/20170409/e9c20f63/attachment.html>