#include <iostream>
#include <iomanip>
#include "vector.h"
#include "chebyshev.h"
#include "flowfield.h"
#include "nsintegrator.h"

int main() {

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

  // Define integration parameters
  const Real dtmax = 0.15;
  const Real dtmin = 0.05;
  const Real CFLmax = 0.90;
  const Real CFLmin = 0.5;
  const Real dT  = 1.0;    // plot interval
  const Real T0  = 100.0;   // start time
  const Real T1  = 200.0;   // stop time

  // Load velocity, modified pressure, and base flow from disk.
  FlowField u("u100");
  FlowField q("q100");
  ChebyCoeff U("U");

  // Get y-domain information from velocity field.
  Real a = u.a();
  Real b = u.b();
  int Ny = u.Ny();

  // Construct Navier-Stoke Integrator
  DNSFlags flags; 
  flags.timestepping = RK3;            // use 3rd-order Runge-Kutta method
  flags.constraint = PressureGradient; // enforce constant pressure gradient

  TimeStep dt((dtmax+dtmin)/2, dtmin, dtmax, dT, CFLmin, CFLmax);
  NSIntegrator dns(u, U, nu, dt, flags, T0);

  ChebyCoeff u00mean(Ny,a,b,Spectral);
  Real dragmean = 0.0;
  int count = 0;

  for (Real t=T0; t<T1; t += dT) {

    // Get kx=kz=0 Fourier component u00(y) and compute drag
    ChebyCoeff u00 = Re(u.profile(0,0,0)); 
    ChebyCoeff du00dy = diff(u00);
    Real drag = nu*du00dy.eval_a();
    
    u00mean += u00;
    dragmean += drag;
    ++count;

    // Save stuff
    string time = i2s(int(t));
    u00.save("uprofile00_"+time);
    Re(u.profile(1,2,0)).save("uprofile12_"+time); 

    cout << "===============================================" << endl;
    cout << "         t == " << t << endl;
    cout << "        dt == " << dt << endl;
    cout << "       CFL == " << dns.CFL() << endl;
    cout << "L2Norm2(u) == " << L2Norm2(u) << endl;
    cout << "      drag == " << drag << endl;

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

    if (dt.adjust(dns.CFL())) {
      cout << "adjusting timestep" << endl;
      dns.reset(nu, dt);
    }
  }

  // Compute means
  dragmean /= count;
  u00mean /= count;
  
  // Fourier-transform u00mean
  u00mean.makePhysical();
  u00mean.save("u00mean");
  cout << "mean drag == " << dragmean << endl;
}