sandbox/easystab/stretching_formula.m

    This code is a validation of the stretching formula for free_surface_mapping.m. In this code, we do not use the numerical mapping map2D.m: we want the analytical formula because the stretching is itself one of the unknown of the system.

    clear all; clf; format compact
disp('%%%%%%%%%')

% parameters
Lx=1;
Ly=1;
Nx=20;
Ny=20;
method='cheb';

% differentiation
[d.x,d.xx,d.wx,x]=dif1D(method,0,Lx,Nx,5);
[d.y,d.yy,d.wy,y]=dif1D(method,0,Ly,Ny,5);
[D,l,X,Y,Z,I,NN]=dif2D(d,x,y);

    Stretching function

    We introduce a stretching of the mesh \displaystyle \bar{x}=x, \bar{y}=\eta(x)y where x and y are the rectangular computational coordinates, and \bar{x} and \bar{y} are the physical coordinates. In this code, everything that is related to the physical domain has subscript p.

    Also, in the doce, we have te variable e which is a vector with Nx elements, and E which is an array with N_y\timesN_x elements, where we copy e for every grid point of the mesh.

    % stretching function eta
e=1-0.3*sin(x*pi);
ex=d.x*e; exx=d.xx*e;
E=ones(Ny,1)*e';E=E(:);
Ex=ones(Ny,1)*ex'; Ex=Ex(:);
Exx=ones(Ny,1)*exx'; Exx=Exx(:);

% physical mesh
Xp=X;
Yp=Y.*reshape(E,Ny,Nx);

    The mapping

    The first useful thing is to get the expression of the derivative of the computational coordinates with respect to the physical coordinates \displaystyle \begin{array}{l} x_\bar{x}=1\\ y_\bar{x}=-y\eta^{-1}\eta_x\\ x_\bar{y}=0\\ y_\bar{y}=\eta^{-1}\\ \end{array}

    Then we use the composition of derivatives to get the expression of the differentiation in physical space in terms of the computational differentiation matrices \displaystyle \begin{array}{l} \phi(x,y)_{\bar{x}}=\phi_x x_\bar{x}+\phi_y y_\bar{x}=\eta^{-1}\phi_y\\ \phi(x,y)_{\bar{y}}=\phi_x x_\bar{y}+\phi_y y_\bar{y}=\phi_x-y\eta^{-1}\eta_x\phi_y\\ \end{array} We do the same operations for the expression of the Laplacian \displaystyle \phi_{\bar{x}\bar{x}}+\phi_{\bar{y}\bar{y}}=(\phi_\bar{x})_\bar{x}+(\phi_\bar{y})_\bar{y} and we get \displaystyle \phi_{\bar{x}\bar{x}}+\phi_{\bar{y}\bar{y}}= \begin{array}{l} \phi_{xx}(1)\\ +\phi_{xy}(-2y\eta^{-1}\eta_x)\\ +\phi_{yy}(y^2\eta^{-2}\eta_x^2+\eta^{-2})\\ +\phi_y(y[2\eta^{-2}\eta_x^2-\eta^{-1}\eta_{xx}])) \end{array} 

    % Build physical space differentiation matrices
Dp.x=D.x ...
    -diag(Y(:).*E.^-1.*Ex)*D.y;
Dp.y=diag(E.^-1)*D.y;

Dp.lap=D.xx ...
    +diag(-2*Y(:).*E.^-1.*Ex)*(D.x*D.y) ...
    +diag(Y(:).^2.*E.^-2.*Ex.^2+E.^-2)*D.yy ...
    +diag(Y(:).*(2*E.^-2.*Ex.^2-E.^-1.*Exx))*D.y;

% Validation
phi=cos(pi*Xp).*sin(pi*Yp);
phix=-pi*sin(pi*Xp).*sin(pi*Yp);
phixx=-pi^2*cos(pi*Xp).*sin(pi*Yp);
phiy=pi*cos(pi*Xp).*cos(pi*Yp);
phiyy=-pi^2*cos(pi*Xp).*sin(pi*Yp);

mesh(Xp,Yp,phi); hold on; 
plot3(x,Ly*e,0*x,'k-');
view(112,38); hold off
xlabel('x'); ylabel('y'); zlabel('phi'); title('phi');
legend('phi','e')

err_x=norm(phix(:)-Dp.x*phi(:))
err_y=norm(phiy(:)-Dp.y*phi(:))
err_lap=norm(phixx(:)+phiyy(:)-Dp.lap*phi(:))

set(gcf,'paperpositionmode','auto');
print('-dpng','-r80','stretching_formula.png');

    Which gives the reassuring screen output

    err_x =
   5.1937e-08
err_y =
   3.8849e-13
err_lap =
   1.8580e-05

    and the figure Figure