function couette

% Define gridsize
Nx=32;
Ny=33;
Nz=32;

% Define box size
Lx=1.75*pi;
a= -1.0;
b=  1.0;
Lz=1.2*pi;

% Define flow parameters
Reynolds = 400.0;
nu = 1.0/Reynolds;
dPdx = 0.0;
Ubulk = 0.0;

% Define integration parameters
n = 25;      % take fifty steps between printouts
dt = 0.04;   % integration timestep
T  = 10.0;  % integrate from t=0 to t=T

% Define DNS parameters 
flags=DNSFlags(...
					'PressureGradient',... %constraint
					'SBDF3',... 			%timestepping 
					'CNRK2',... 			%initstepping
					'Rotational',... 	%nonlinearity
					'DealiasXZ',... 	%dealiasing
					'true',... 			%taucorrection
					'PrintTicks',...	%verbosity
					dPdx,...  			%dPdx
					Ubulk... 				%Ubulk
					);

% Define size and smoothness of initial disturbance
spectralDecay = 0.3;   
magnitude  = 0.01;   
kxmax = 3;
kzmax = 3;

% Construct base flow for plane Couette: U(y) = y
U=chebypoints(Ny,a,b,'Physical'); %U is a ChebyCoeff object
 
% Construct data fields: 3d velocity and 1d pressure
printf('building velocity and pressure fields...');
u=FlowField([Nx,Ny,Nz],3,Lx,Lz,a,b);
q=FlowField([Nx,Ny,Nz],1,Lx,Lz,a,b);
printf('done.');

% Perturb velocity field
addPerturbations(u,kxmax,kzmax,1.0,spectralDecay);
scaleby(u,magnitude/L2Norm(u));

% Construct Navier-Stoke integrator, set integration method
printf('building DNS...');
dns=DNS(u,U,nu,dt,flags); 
printf('done');

for t=0:n*dt:T
	printf('t = %g',t);
	printf('CFL(dns) = %g',CFL(dns));
	printf('L2Norm(u) = %g',L2Norm(u));
	printf('divNorm(u) = %g',divNorm(u));
	advance(dns,u,q,n);

end;