Warning, /verification/advect_xz/input/gendata.m is written in an unsupported language. File is not indexed.
view on githubraw file Latest commit add29e06 on 2018-01-31 20:35:05 UTC
c2bfce3ff6 Jean*0001 % This is a matlab script that generates the input data
0002 prec='real*8';
0003 ieee='b';
0004 w2file=1; %- do write to file (1=y/0=no)
0005
0006 %- Dimensions of grid
0007 nx=20;
0008 ny=1;
0009 nz=20;
0010
0011 %- Horizontal & vertical resolution (m):
0a0cf9eb41 Jean*0012 dx=1.e4; dy=dx; dz=1.e2;
c2bfce3ff6 Jean*0013 % full size of the domain:
0014 Lx=nx*dx;
0015 H=nz*dz ;
0016
0017 %- grid point coordinate :
0018 xf=[0:nx]*dx;
0019 xc=(xf(1:nx)+xf(2:nx+1))/2;
0020 zc=-dz/2:-dz:-H;
0021 zf=0:-dz:-H;
0022
0023 wf=' write file: ';
0024 %- bathymetry :
0025 zb=-H*ones(nx,1);
0026 zb=-H*(1.-xc'/Lx/2);
0027
0028 if w2file,
0029 bname='bathy_slope.bin';
0030 fid=fopen(bname,'w',ieee); fwrite(fid,zb,prec); fclose(fid);
0031 fprintf([wf,bname,'\n']);
0032 end
0033
0034 %- partial cell filling (hFac):
0035 hFac=ones(nx,1)*zf(1,[2:nz+1]);
0036 hFac=max(hFac,zb*ones(1,nz));
0037 hFac=ones(nx,1)*zf(1,[1:nz])-hFac;
0038 hFac=hFac/dz; hFac=max(0,hFac);
0039
0040 rhFc=reshape(hFac,[nx*nz 1]);
0041 [IK]=find(rhFc > 0); rhFc(IK)=1./rhFc(IK);
0042 rhFc=reshape(rhFc,[nx nz]);
0043
0044 %- horizontal flow field is defined from stream-function psi :
0045 psi=zeros(nx+1,nz+1);
0a0cf9eb41 Jean*0046 psi1=sin(pi*xf'/Lx)*ones(1,nz+1);
c2bfce3ff6 Jean*0047 zu=zeros(nx+1,1); zu(2:nx)=max(zb(2:nx),zb(1:nx-1));
0048 rHu=zeros(nx+1,1); rHu(2:nx)=-1./zu(2:nx);
0049 psi2=max( ones(nx+1,1)*zf, zu*ones(1,nz+1) );
0050 psi2=psi2.*(rHu*ones(1,nz+1));
0051 psi2=sin(pi*psi2);
0a0cf9eb41 Jean*0052 psi=H*psi1.*psi2;
c2bfce3ff6 Jean*0053
0054 uTrans=psi(:,[1:nz])-psi(:,[2:nz+1]);
0055 u0=uTrans(1:nx,:).*rhFc/dz;
0056
0a0cf9eb41 Jean*0057 %- add small, vertically uniform, divergent component:
0058 ud=psi1(:,1).*rHu;
0059 ud=ud*H/160.;
0060 u1=ud(1:nx,1)*ones(1,nz);
0061 u1(find(hFac==0.))=0.;
0062
c2bfce3ff6 Jean*0063 if w2file,
0064 uname='Uvel.bin';
0065 fid=fopen(uname,'w',ieee); fwrite(fid,u0,prec); fclose(fid);
0066 fprintf([wf,uname,'\n']);
0a0cf9eb41 Jean*0067 uname='Udiv.bin';
0068 fid=fopen(uname,'w',ieee); fwrite(fid,u0+u1,prec); fclose(fid);
0069 fprintf([wf,uname,'\n']);
c2bfce3ff6 Jean*0070 end
0071
0072 %- Initial Temperature : Gaussian bump :
0073 t0=zeros(nx,nz);
0074 % position of the center
0075 xM=xf(5); zM=zf(6);
0076 % size of the patch
0077 dX=2*dx; dH=2*dz;
0078 % amplitude:
0079 tA=exp(-0.5*(35/20)^2);
0080 %
0081 r1=(xc'-xM)/dX; r2=(zc-zM)/dH;
0082 r1=r1.*r1; r2=r2.*r2;
0083 rD=r1*ones(1,nz)+ones(nx,1)*r2;
0084 t0=tA*exp(-rD/2);
0085
0086 if w2file,
0087 tname='Tini_G.bin';
0088 fid=fopen(tname,'w',ieee); fwrite(fid,t0,prec); fclose(fid);
0089 fprintf([wf,tname,'\n']);
0090 end
0091
0092 %return
0093
0094 %- make plots to check:
0095 figure(1);clf;
0096 subplot(211)
0a0cf9eb41 Jean*0097 %imagesc(xf,zf,psi1'); set(gca,'YDir','normal'); colorbar;
c2bfce3ff6 Jean*0098 imagesc(xf,zf,psi'); set(gca,'YDir','normal'); colorbar;
0099 hold on ; plot(xc,zb,'b-'); hold off ; grid
0100 title('Stream-Function [m^2/s]');
0101 subplot(212)
0a0cf9eb41 Jean*0102 %imagesc(xf(1:nx),zc,u1'); set(gca,'YDir','normal'); colorbar;
c2bfce3ff6 Jean*0103 imagesc(xf(1:nx),zc,u0'); set(gca,'YDir','normal'); colorbar;
0104 hold on ; plot(xc,zb,'r-'); hold off ; grid
0105 title('U velocity [m/s]');
0106
0107 figure(2);clf;
0108 subplot(211)
0109 imagesc(xc,zc,rhFc'); set(gca,'YDir','normal'); colorbar;
0110 hold on ; plot(xc,zb,'k-'); hold off ; grid
0111 title('Recip hFac');
0112 subplot(212)
0113 imagesc(xc,zc,t0'); set(gca,'YDir','normal'); colorbar;
0114 hold on ; plot(xc,zb,'r-'); hold off ; grid
0115 title('Initial Theta');
0116
0117 return
