[Nek5000-users] intp_rstd; coarse to fine and back

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

Hi Sandeep,


intp_rstd(....,1)


does not interpolate back.


It takes the transpose of the forward interpolation operator.


WHY?


Because if I want to integrate two functions using quadrature on

a fine mesh, I need a routine that interpolates _each_ function to the

fine mesh.


Thus, to evaluate


     (v,u) := \int  v u dx   =  v^T B u


I would use a quadrature rule on M gauss points, where M > N (typically).


If u = [ u0 u1 ... uN ]^T is the set of basis coefficients for u, then evaluating


u(x) on points xi_k, with xi_k, k=1,...,M being M gauss points, requires:


      u(xi_k) = \sum_j=0^N    h_j (xi_k) u_j


Let uf = [ u(xi_1)  u(xi_2) ... u(xi_M) ]^T be the vector of u values on the

fine points.


Then:


       uf = J u


where J_{kj} :=  h_j(xi_k)


and h_j(x) is the jth cardinal Lagrange polynomial of degree N on the

spectral element nodal points.


Let  rho_k, k=1,...,M  be the quadrature weights for the M-point gauss rule.


Then,


(v,u) ==   \sum_k  vf_k  rho_k  uf_k


          =    (Jv)^T  D  Ju


          =    v^T  J^T D  J u


          =    v^T     B        u


where B:=J^T D J, and D = diag (rho_k), k=1,...,M.


So, when evaluating a matrix-vector product of the form


            w = B u


we could (and often do) express this as:


             w1  =    Ju


              w2  =   D w1


              w  = J^T  w1



Now, done.


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


hth,


Paul




________________________________
From: Nek5000-users <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, June 4, 2018 9:24:20 AM
To: nek5000-users at lists.mcs.anl.gov
Subject: [Nek5000-users] intp_rstd; coarse to fine and back

Hello,
To know more about the subroutines used for dealiasing, I perform a simple exercise.
I take a simple sine function and project it onto the fine mesh and back to the coarse using intp_rstd.
When I compare the values of the processed array with the original array, it is scaled by a factor of 1.25204, ie, I get a function like 1.25204*sin(2*x)
Could you please tell me which factor am I missing? Which jacobian has to be multiplied?
Attached is the usr file and below are the operations I perform in short.

      sinax(i,1,1,1)      = sin(2.*x)

     call intp_rstd(uf1,u1,nx1,nxd,if3d,0) ! 0 --> forward
      nxyzd = lxd*lyd*lzd
         do i=1,nxyzd
            tr1   = uf1(i)
            wf(i) = tr1*rx(i,1,e)+tr1*rx(i,2,e)+tr1*rx(i,3,e)
         enddo

      call intp_rstd(w1,wf,nx1,nxd,if3d,1) ! 1 --> back to coarse

      nxyz = lx1*ly1*lz1           !
         do i=1,nxyz               ! to physical space.
            bi = 1./((bm1(i,1,1,e)))   !
            w1(i) = bi*w1(i)
            jacm=jacmi(i,e) ! I dont know whether to use jacmi(i,e) or not
         enddo
Thank you in advance.

Sandeep


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