%{ Validation: Comparing the expression of the analytical Jacobian and the numerical Jacobian for the functions on a mapped domain. We are only interested in the part of the Jacobian where the expression of the stretching plays a role. This is a follow-up of the code [stretching_formula.m]() and to be used in [domain_derivative_mapping.m](). %} disp('%%%%%%%') clear all; clf; % parameters Nx=10; % number of grid points in x Ny=10; % number of grid points in y Lx=1; % the length in x of the computational domain Ly=1; % the length in x of the computational domain eps=1e-8; % Small perturbation for computing the numerical Jacobian method='cheb'; % discretization method pts=5; % number of points in the stencil in case of finite differences % differentiation [d.x,d.xx,d.wx,x]=dif1D(method,0,Lx,Nx,pts); [d.y,d.yy,d.wy,y]=dif1D(method,0,Ly,Ny,pts); [D,l,X,Y,Z,I,NN]=dif2D(d,x,y); % An operator to transform from the 1D x grid to the 2D grid T=kron(speye(Nx,Nx),ones(Ny,1)); yy=Y(:); l.phi=(1:NN)'; l.e=NN+[1:Nx]'; n=NN+Nx; % base state e=1+0.2*sin(pi*x/Lx); phi=cos(pi*X).*cos(pi*Y); q=[phi(:); e(:)]; % Computing Jacobian with respect to e A=zeros(3*NN,Nx); for gre=0:Nx % Perturbation of one degree of freedom if gre>0; qq=q+((1:n)==(l.e(gre)))'*eps; else qq=q; end % computational domain derivatives phi=qq(l.phi); phiy=D.y*phi; phixx=D.xx*phi; phiyy=D.yy*phi; phixy=D.x*phiy; e=qq(l.e); ex=d.x*e; exx=d.xx*e; E=T*e; Ex=T*ex; Exx=T*exx; %{ # The physical space differentiation matrices Here instead of using the general function [map2D.m]() to get the differentiation matrices on the deformed domain, we build them from the analytical mapping as derived and tested in [stretching_formula.m#the-mapping](). %} % Physical space diferentiation matrices Dp.lap=D.xx ... +spd(-2*yy.*E.^-1.*Ex)*(D.x*D.y) ... +spd(yy.^2.*E.^-2.*Ex.^2+E.^-2)*D.yy ... +spd(yy.*(2*E.^-2.*Ex.^2-E.^-1.*Exx))*D.y; Dp.x=D.x+spd(-yy.*E.^-1.*Ex)*D.y; Dp.y=spd(E.^-1)*D.y; % Compute the physical space derivatives phiplap=Dp.lap*phi; phipx=Dp.x*phi; phipy=Dp.y*phi; % nonlinear function feps=[phiplap; phipx; phipy]; % store Jacobian if gre==0; % The analytical expressions for the Jacobian f=feps; %{ # The analytical expression of the Jacobian When we just had to get the Jacobian of derivatives it was easy because these are linear operators, thus the Jacobian is the same as the operator itself. Now it is a bit more complicated because the expression of the derivative depends on the domain stretching $\eta$, and thus the derivatives are nonlinear functions of $\eta$^[cf jerome's notebook c107p130-131]. In the following we denote $$\partial_\bar{x}^{\eta}$$ the physical-space differentiation operator based on a domain stretched by $\eta$. And we are interested in how this thing depends on $\eta$ for a small perturbation $\tilde{\eta}$. # First derivative in $y$ We have $$\partial^\eta_\bar{y}\phi=\eta^{-1}\phi_y$$ thus considering the perturbations of $\phi$ and $\eta$ we have $$\partial^{\eta+\tilde{\eta}}_\bar{y}(\phi+\tilde{\phi}) =(\eta+\tilde{\eta})^{-1}(\phi_y+\tilde{\phi}_y)\approx \eta^{-1}\phi_y+\eta^{-1}\partial_y\tilde{\phi}-\eta^{-2}\phi_y\tilde{\eta}$$ which we can write in matrix form as $$\partial^{\eta+\tilde{\eta}}_\bar{y}(\phi+\tilde{\phi}) \approx \partial^\eta_\bar{y}\phi+ \begin{pmatrix} \partial^\eta_\bar{y} &-\eta^{-2}\phi_y\end{pmatrix} \begin{pmatrix} \tilde{\phi}\\ \tilde{\eta} \end{pmatrix}$$ # First derivative in $x$ We have $$\partial^\eta_\bar{x}\phi=\phi_x-y\eta^{-1}\eta_x\phi_y$$ then considering the perturbations of $\phi$ and $\eta$ we have $$\partial^{\eta+\tilde{\eta}}_\bar{x}(\phi+\tilde{\phi}) \approx \partial^{\eta}_\bar{x}\phi+ \begin{pmatrix} \partial^\eta_\bar{x} &-y\phi_y(\eta^{-1}\partial_x-\eta^{-2}\eta_x) \end{pmatrix} \begin{pmatrix} \tilde{\phi}\\ \tilde{\eta} \end{pmatrix}$$ # For the laplacian Now for the physical-space Laplacian $$\bar{\Delta}^{\eta}=\partial^\eta_{\bar{xx}}\partial^\eta_{\bar{yy}}$$ we have $$\bar{\Delta}^{\eta+\tilde{\eta}}(\phi+\tilde{\phi}) =\bar{\Delta}_\eta \phi+ \begin{pmatrix} \bar{\Delta}^{\eta} & \phi_{xx}A+\phi_{yy}B+\phi_{xx}C \end{pmatrix} \begin{pmatrix} \tilde{\phi}\\ \tilde{\eta} \end{pmatrix}$$ with $A, B$ and $C$ the Jacobians of the functions of $\eta$ $$\begin{array}{l} a(\eta)=-2y\eta^{-1}\eta_x\\ b(\eta)=y^2\eta^{-2}\eta_x^2+\eta^{-2}\\ c(\eta)=y(2\eta^{-2}\eta_x^2-\eta^{-1}\eta_{xx} \end{array}$$ After some algebra, we find $$\begin{array}{l} A=-2y(\eta^{-1}\partial_x-\eta^{-2}\eta_x)\\ B=2y^2\eta^{-2}(\eta_x\partial_x-\eta^{-1}\eta_x^2)-2\eta^{-3}\\ C=y[ 4\eta^{-2} (\eta_x\partial_x-\eta^{-1}\eta^2_{x})-\eta^{-1}(\partial_{xx}-\eta^{-1}\eta_{xx})] \end{array}$$ Now the difficulty with these expressions is that $y$ is defined for everypoint of the grid, but $\eta$ is only define on a 1D grid, so we need to use an operator that multiplies the size *Nx* vector of $\eta$ to give the grid equivalent of $\eta$. This is the operator *T* built above. %} AA=-2*spd(yy)*T*(diag(e.^-1)*d.x-diag(e.^-2.*ex)); BB=2*(spd(yy.^2)*T*(diag(e.^-2.*ex)*d.x-diag(e.^-3.*ex.^2))-T*diag(e.^-3)); CC=spd(yy)*T*(diag(4*e.^-2.*ex)*d.x-diag(e.^-1)*d.xx+diag(-4*e.^-3.*ex.^2+e.^-2.*exx)); Dp.lapphie=spd(phixy)*AA+spd(phiyy)*BB+spd(phiy)*CC; Dp.xphie=spd(yy.*phiy)*T*(diag(-e.^-1)*d.x+diag(e.^-2.*ex)); Dp.yphie=spd(phiy)*T*diag(-e.^-2); else % One column of the numerical Jacobian A(:,gre)=(feps-f)/eps; end end % Validation: compare numerical and analytical Jacobians err_lap=norm(A(1:NN,:)-Dp.lapphie) err_x=norm(A(NN+(1:NN),:)-Dp.xphie) err_y=norm(A(2*NN+(1:NN),:)-Dp.yphie) % Show the function phi and e Yp=reshape(E.*Y(:),Ny,Nx); mesh(X,Yp,reshape(phi,Ny,Nx)); hold on plot(x,e,'k-',[0 Lx],[Ly Ly],'k--'); view(103,66); xlabel('x');ylabel('y');zlabel('h'); set(gcf,'paperpositionmode','auto') print('-dpng','-r80','jacobian_formula.png') %{ This gives the reassuring screen output err_lap = 9.5728e-04 err_x = 9.9049e-06 err_y = 1.0995e-06 # The figure The figure just shows the reference state $\phi$ and $\eta$ used to test the analytical expressions of the Jacobian. ![](jacobian_formula.png) # Exercices/Contributions * Please * Please %}