<html dir="ltr">
<head>
<meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
<style type="text/css" style="">
<!--
p
{margin-top:0;
margin-bottom:0}
-->
</style><style type="text/css" id="owaParaStyle"></style>
</head>
<body dir="ltr" fpstyle="1" ocsi="0">
<div style="direction: ltr;font-family: Tahoma;color: #000000;font-size: 10pt;">Dear Juan,
<div><br>
</div>
<div>The usage model for intp_rstd() is to map back and forth between a std grid and a fine grid,</div>
<div>in the context of Galerkin projection. Here is the basic idea:</div>
<div><br>
</div>
<div>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</div>
<div>is the closest element in PN to u*, then you have to choose a norm. </div>
<div><br>
</div>
<div>Let's suppose you are choosing the L2 norm. Let's also assume u*=vx^2, which is a polynomial of</div>
<div>degree 2N.</div>
<div><br>
</div>
<div>If u is the closest element of PN to u*, it means that (u*-u) \perp PN, which is to say,</div>
<div><br>
</div>
<div> (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].</div>
<div><br>
</div>
<div>Assuming there are n basis functions in PN, we have:</div>
<div><br>
</div>
<div> u = sum_j=1^n u_j phi_j (X)</div>
<div><br>
</div>
<div>and thus we have n unknown basis coefficients, u_j, j=1,...n.</div>
<div><br>
</div>
<div>Moreover, the projection (i.e., best-fit) statement gives us n equations:</div>
<div><br>
</div>
<div> (phi_i , u-u* ) = 0 , i=1,...,n</div>
<div><br>
</div>
<div>or,</div>
<div><br>
</div>
<div> (phi_i, u ) = (phi_i , u*) , i=1,...,n</div>
<div><br>
</div>
<div>Bescause the inner product<span style="font-size: 10pt;"> (f,g) := \int_Omega f g dX </span><span style="font-size: 10pt;">is linear, we </span><span style="font-size: 10pt;"> can right the term on the left as</span></div>
<div><span style="font-size: 10pt;"><br>
</span></div>
<div><span style="font-size: 10pt;"> sum_j (phi_i, phi_j ) u = (phi_i , u*) =: b_i , i=1,...,n</span></div>
<div><span style="font-size: 10pt;"><br>
</span></div>
<div>so B u = b</div>
<div><br>
</div>
<div>where B_{ij} := (phi_,phi_j) := \int_Omega phi_i , phi_j dX.</div>
<div><br>
</div>
<div>Thus, we have two types of integrals to compute, B_{ij} and b_i. Let's take care of b_i first:</div>
<div><br>
</div>
<div> b_i = \int_Omega v u* dX</div>
<div><br>
</div>
<div>If u* \in P_2N and v \in P_N, then we need a quadrature rule that is exact for all polynomials of</div>
<div>degree M \leq 3N. Gauss (not Gauss Lobatto) quadrature on (3M/2+1) points is exact for all</div>
<div>polynomials of degree 2(3M/2+1) - 1 = 3M+1, so that will suffice. To compute b_i, we thus</div>
<div>do:</div>
<div><br>
</div>
<div> b_i = sum_k=1^m phi_i (\xi_k) rho_k vx(\xi_k) vx(\xi_k)</div>
<div><br>
</div>
<div>where \xi_k are the Gauss quadrature points in the domain, m=number of Gauss points,</div>
<div>rho_k are the Mth-order Gauss quadrature weights.</div>
<div><br>
</div>
<div>Note that, it is important to interpolate vx to the \xi_k, and not u*, because information is</div>
<div>otherwise already lost. Note that if</div>
<div><br>
</div>
<div> vx(X) = sum_j \phi_j(X) vx_j</div>
<div><br>
</div>
<div>then</div>
<div><br>
</div>
<div> vx(\xi_k) = \sum \phi_j(\xi_k) vx_j = J vx</div>
<div><br>
</div>
<div>where J_{kj} := \phi_j(\xi_k) is the interpolation matrix from the Nth-order Gauss Lobatto (GL)</div>
<div>points to the Mth-order Gauss (G) quadrature points. Careful inspection of of the above expression</div>
<div>shows that the rhs ( b = {b_i} ) is given by:</div>
<div><br>
</div>
<div> b = J^T B_M V (J vx)</div>
<div><br>
</div>
<div>where V = diag ( vx(\xi_k) ) is just a diagonal matrix with vx on the fine quadrature mesh.</div>
<div><br>
</div>
<div>So, ignoring momentarily the diagonal matrix V, the process of generating the right hand side,</div>
<div>with entries on the Nth-order Gauss points is, in words:</div>
<div><br>
</div>
<div> Interpolate vx onto fine mesh: vx --> J vx</div>
<div><br>
</div>
<div> Collocate with V and the diagonal matrix of quadrature weights, B_M</div>
<div><br>
</div>
<div> Map back to the Nth-order GL mesh: J^T ( B_M V (J vx ) )</div>
<div><br>
</div>
<div>OK. So,</div>
<div><br>
</div>
<div> intp_rstd(...,0) applies J.</div>
<div><br>
</div>
<div> intp_rstd(...,1) applies J^T</div>
<div><br>
</div>
<div>Now, you still haven't finished the projection process.</div>
<div><br>
</div>
<div>To get u, you must do the following:</div>
<div><br>
</div>
<div> call col3(u,binvm1,rhs,n)</div>
<div><br>
</div>
<div>which applies B^{-1} to get u (from, as it were, u*= vx * vx ).</div>
<div><br>
</div>
<div><br>
</div>
<div>As a final finishing touch, the "col3" statement will be exact if the product u*v is in P_2N-1. Unfortunately,</div>
<div>it is in P_2N, because u and v are both polynomials of degree N.</div>
<div><br>
</div>
<div>You have the option of computing the non diagonal matrix and solving that system via preconditioned</div>
<div>conjugate gradients, or of just allowing that small discrepancy. If you've come this far, I would </div>
<div>recommend the former. I have code that can do that for you.</div>
<div><br>
</div>
<div>Paul</div>
<div><br>
</div>
<div><br>
</div>
<div><br>
</div>
<div><br>
</div>
<div><br>
</div>
<div><br>
</div>
<div>
<div style="font-family: Times New Roman; color: #000000; font-size: 16px">
<hr tabindex="-1">
<div id="divRpF73600" style="direction: ltr;"><font face="Tahoma" size="2" color="#000000"><b>From:</b> nek5000-users-bounces@lists.mcs.anl.gov [nek5000-users-bounces@lists.mcs.anl.gov] on behalf of nek5000-users@lists.mcs.anl.gov [nek5000-users@lists.mcs.anl.gov]<br>
<b>Sent:</b> Sunday, April 09, 2017 2:12 AM<br>
<b>To:</b> nek5000-users@lists.mcs.anl.gov<br>
<b>Subject:</b> [Nek5000-users] Problem with intp_rstd<br>
</font><br>
</div>
<div></div>
<div>
<div id="divtagdefaultwrapper" dir="ltr" style="font-size:12pt; color:#000000; font-family:Calibri,Arial,Helvetica,sans-serif">
<p>Dear all,</p>
<p><br>
</p>
<p>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. </p>
<p><br>
</p>
<p></p>
<div> do e=1,nelv ! mapping to fine grid</div>
<div> call intp_rstd(vxd(1,1,1,e),vx(1,1,1,e),nx1,nxd,if3d,0) ! 0 --> forward</div>
<div> enddo</div>
<div><br>
</div>
<div> call vsq(vxd,ntotd) ! calculating vxd**2</div>
<div></div>
<div><br>
</div>
<div> do e=1,nelv ! mapping back to original grid</div>
<div> call intp_rstd(vx(1,1,1,e),vxd(1,1,1,e),nx1,nxd,if3d,1) ! 0 --> backwards</div>
<div> enddo</div>
<br>
<p></p>
<p><span style="font-family:Calibri,Arial,Helvetica,sans-serif,EmojiFont,"Apple Color Emoji","Segoe UI Emoji",NotoColorEmoji,"Segoe UI Symbol","Android Emoji",EmojiSymbols; font-size:16px">Am I doing something wrong?</span><br>
</p>
<p><br>
</p>
<p>Thanks,</p>
<p><br>
</p>
<p>Juan Diego</p>
</div>
</div>
</div>
</div>
</div>
</body>
</html>