#include <iostream>
#include <iomanip>
#include "channelflow/dns.h"
#include "channelflow/flowfield.h"

using namespace std;

int main() {

  // Define gridsize
  const int Nx=32;
  const int Ny=33;
  const int Nz=32;

  // Define box size
  const Real Lx=1.75*pi;
  const Real a= -1.0;
  const Real b=  1.0;
  const Real Lz=1.2*pi;

  // Define flow parameters
  const Real Reynolds = 400.0;
  const Real nu = 1.0/Reynolds;

  // Define integration parameters
  const int n = 25;      // take fifty steps between printouts
  const Real dt = 0.04;  // integration timestep
  const Real T  = 100.0; // integrate from t=0 to t=T

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

  // Construct base flow for plane Couette: U(y) = y
  ChebyCoeff U(Ny,a,b,Physical);
  Vector y = chebypoints(Ny, a,b);
  for (int ny=0; ny<Ny; ++ny) 
    U[ny] = y[ny];
  U.save("U");
  y.save("y");
 
  // Construct data fields: 3d velocity and 1d pressure
  cout << "building velocity and pressure fields..." << flush;
  FlowField u(Nx,Ny,Nz,3,Lx,Lz,a,b);
  FlowField q(Nx,Ny,Nz,1,Lx,Lz,a,b);
  cout << "done" << endl;

  // Perturb velocity field
  u.addPerturbations(kxmax,kzmax,1.0,spectralDecay);
  u *= magnitude/L2Norm(u);

  // Construct Navier-Stoke integrator, set integration method
  cout << "building DNS..." << flush;
  DNSFlags flags; 
  flags.timestepping = SBDF3;
  flags.initstepping = CNRK2;
  flags.nonlinearity = Rotational;
  flags.dealiasing = DealiasXZ;
  flags.constraint  = PressureGradient; // enforce constant pressure gradient

  DNS dns(u, U, nu, dt, flags);
  cout << "done" << endl;

  for (Real t=0; t<=T; t += n*dt) {
    cout << "         t == " << t << endl;
    cout << "       CFL == " << dns.CFL() << endl;
    cout << " L2Norm(u) == " << L2Norm(u) << endl;
    cout << "divNorm(u) == " << divNorm(u) << endl;

    // Write velocity and modified pressure fields to disk
    u.binarySave("u"+i2s(int(t)));
    q.binarySave("q"+i2s(int(t)));

    // Take n steps of length dt
    dns.advance(u, q, n);
    cout << endl;
  }
}