/*==============================================================================

                              ****************************
                              *         P  N  S  2       *
                              ****************************

     A FINITE ELEMENT CODE FOR DYNAMIC INCOMPRESSIBLE FLUID FLOW ANALYSIS
                                  IN 2-D GEOMETRIES

  ==============================================================================

   Copyright (C) 1998 - 2002 Rachid Touzani

   This program is free software; you can redistribute it and/or modify it under
   the terms of the GNU General Public License as published by the Free 
   Software Foundation; Version 2 of the License.

   This program is distributed in the hope that it will be useful, but WITHOUT
   ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
   FOR A PARTICULAR PURPOSE. See the GNU General Public License for more 
   details.

   You should have received a copy of the GNU General Public License 
   along with this program; if not, write to the :

   Free Software Foundation
   Inc., 59 Temple Place - Suite 330
   Boston, MA  02111-1307, USA

  ==============================================================================*/

#include "main.h"

int main(int argc, char *argv[])
{
   int      nb_it, step, m_it=0;
   ofstream lf;
   double   tm = 0;
   double   cfl = 0;
   FDF      *vff=NULL, *pff=NULL, *sff=NULL, *wff=NULL;
   VDF      *vdf=NULL;
   Timer    cpu_time[4];

   if (argc < 2) {
     cout << "\n\nUsage:  pns2  <parameter_file>" << endl;
     return 0;
   }

   IPF data("pns2 - 1.0",argv[1]);
   int verbose = data.Verbose();
   int output_flag = data.Output();
   int save_flag = data.Save();
   int plot_flag = data.Plot();
   int init_flag = data.Init();
   double deltat = data.TimeStep();
   double max_time = data.MaxTime();
   double dens = data.DoublePar(1);
   double visc = data.DoublePar(2);
   double tol = data.Tolerance();

   lf.open(data.AuxFile(1),ios::out);
   lf << "=================================================================================\n\n";
   lf << "                                   PNS2\n\n";
   lf << "                       A Finite Element Code for\n";
   lf << "             2-D Dynamic Incompressible Fluid Flow Analysis\n\n";
   lf << "                          V E R S I O N   2.1\n\n";
   lf << "=================================================================================\n\n";

//----------
// Read data
//----------


// Read Mesh data
   if (verbose > 1)
     cout << "Reading mesh data ...\n";
   Mesh ms(data.MeshFile(1));
//   ms.NumberEquations(1);
   int nb_dof = ms.Dim();
   int nb_eq = ms.NbNodes();

// Allocate memory for matrices and vectors
   if (verbose > 1)
     cout << "Allocating memory for matrices and vectors ...\n";
   NodeVect<> uf(ms,nb_dof,"Velocity");
   NodeVect<> pf(ms,1,"Pressure");
   NodeVect<> psi(ms,1,"Stream");
   Vect<> u0(2*nb_eq), u(2*nb_eq), v(2*nb_eq), c(2*nb_eq);
   Vect<> ub(2*ms.NbElements()), p(nb_eq), q(nb_eq);

   cout << "Running pns2 ...\n\n";
   cout << "Read execution information in file : ";
   cout << data.AuxFile(1) << endl;
//   if (data.Data())
     vdf = new VDF(ms,data.DataFile());

// Print Mesh data
   if (verbose > 5)
     cout << ms;

// Read initial condition, boundary conditions, body and boundary forces

   int nb_step = int((max_time-tm)/deltat);
   if (save_flag > nb_step || save_flag == 0)
     save_flag = nb_step;
   if (plot_flag > nb_step)
     plot_flag = nb_step;
   if (plot_flag) {
     vff = new FDF(data.AuxFile(2),FDF_WRITE);
     pff = new FDF(data.AuxFile(3),FDF_WRITE);
     sff = new FDF(data.AuxFile(4),FDF_WRITE);
     wff = new FDF(data.AuxFile(5),FDF_WRITE);
   }


// Print File Information
// ----------------------

   lf << "FILE INFORMATION :\n\n";
   lf << "Mesh file :                            " << data.MeshFile(1) << endl;
   if (data.Init()==GET_FROM_FDF)
     lf << "Initial Data file :                    " << data.InitFile() << endl;
   if (data.BF()==GET_FROM_FDF)
     lf << "Body Force file :                      " << data.BFFile() << endl;
   if (data.SF()==GET_FROM_FDF)
     lf << "Boundary Traction file :               " << data.SFFile() << endl;
   if (plot_flag) {
     lf << "Velocity saved in file :               " << data.AuxFile(2) << endl;
     lf << "Pressure saved in file :               " << data.AuxFile(3) << endl;
     lf << "Stream Function saved in file :        " << data.AuxFile(4) << endl;
     lf << "Vorticity saved in file :              " << data.AuxFile(5) << endl;
     lf << "Intermedite Solutions saved in file  : " << data.AuxFile(6) << endl;
   }
   lf << endl << endl << "Viscosity       : " << visc << endl;
   lf << "Density         : " << dens << endl;
   lf << "Time Step       : " << deltat << endl << endl;


// Calculate Pressure and Velocity Matrices
// ----------------------------------------

   cpu_time[0].Start();
   Vect<> M(nb_eq);
   if (verbose > 1)
     cout << "Calculating Pressure Matrix ...\n";
   SpMatrix<> Ap(ms,NODE_DOF,1);
   PMatrix(deltat/dens,ms,Ap,M);
   ILUPrec< double,SpMatrix<> > *Pp = new ILUPrec< double,SpMatrix<> >(Ap);
   cpu_time[0].Stop();
   SpMatrix<> Av(ms);
   VMatrix(deltat,dens,visc,ms,Av);
   ILUPrec< double,SpMatrix<> > Pv(Av);

// ====================
// Loop over time steps
// ====================

   u0 = u;
   if (data.Init()==GET_FROM_VDF)
     vdf->Get(INITIAL,u);

   InitialV(deltat/dens, ms, Ap, M, *Pp, v, u, ub, p, q);
   lf.setf(ios::scientific);
   lf << "\n\nTIME INTEGRATION HISTORY :\n";
   lf << "\nTime Step      Time         Velocity Norm      cfl            Stationary Conv.\n\n";

   for (step=1; step<=nb_step; step++) {

      tm += deltat;
      if (verbose > 1)
        cout << "\n---------------------------------------------------------------------" << endl;
      cout << "Time step " << step << " ... Time = " << tm << endl;


//    Solve Momentum Equations
//    ------------------------

      if (verbose > 1)
        cout << "Solve Momentum Equations ..." << endl;
      cpu_time[1].Start();
      Momentum(step, tm, deltat, ms, *vdf, u, v, ub, c, Av, Pv, p, dens, visc, cfl, verbose);
      if (verbose > 0)
        cout << "CFL = " << cfl << endl;
      cpu_time[1].Stop();

//    Solve Pressure Equation and Update Velocity
//    -------------------------------------------

      if (verbose > 1)
        cout << "Calculating Pressure ..." << endl;
      cpu_time[2].Start();
      nb_it = CalPres(deltat, ms, Ap, *Pp, u, ub, p, q);

      if (verbose > 0)
        cout << "Nb. of CG Iterations for pressure   : " << nb_it << endl;
      cpu_time[2].Stop();

      if (verbose > 1)
        cout << "Updating Velocity ..." << endl;
      cpu_time[3].Start();
      UpdateV(ms, deltat/dens, M, q, u, ub);
      cpu_time[3].Stop();

//    Print Solution Information
//    --------------------------

      m_it += nb_it;
      u0 -= u;
      double diff = u0.NormMax()/deltat;
      double uu = u.NormMax();
      lf << setw(6) << step << "      " << setprecision(2) << setw(10) << tm;
      lf << "      " << setprecision(4) << setw(8) << uu << "      ";
      lf << setprecision(4) << setw(8) << cfl << "   ";
      lf << "      " << setprecision(6) << setw(10) << diff << endl;
      if (verbose > 1) 
        cout << "Velocity Max-Norm    : " << uu << endl;
      if (verbose > 0)
        cout << "Discrepancy          : " << diff << endl;
      u0 = u;
      if (diff < tol) {
        lf << "\n\n*** Convergence to Stationary Solution ***\n";
        cout << "\n---------------------------------------------------------------------\n";
        cout << "Convergence to Stationary Solution." << endl;
        step = nb_step;
      }

//    Output and/or Store solution

      uf.Fill(u); uf.Time(tm); psi.Time(tm);
      pf.Fill(p,1); pf.Time(tm);

      if (output_flag > 0 && output_flag%step==0) {
        cout << uf;
        cout << p;
      }

//    Save in plot file

      if (plot_flag > 0 && step%plot_flag==0) {
        vff->Put(uf);
        pff->Put(pf);
        StreamF(ms,u,psi);
        psi.Name("Stream");
        sff->Put(psi);
        Vorticity(dens/deltat,ms,u,M,psi);
        psi.Name("Vorticity");
        wff->Put(psi);
      }

//    Save in restart file

      if (save_flag > 0 && save_flag%step==0) {
        FDF ffz(data.AuxFile(6),FDF_WRITE);
        ffz.Put(u,tm);
        ffz.Put(ub,tm);
        ffz.Put(p,tm);
      }
   }

   lf << "\n\nMean Nb. of Iterations per Time Step : " << m_it/step << endl;
   lf << "\n\nCPU TIMES (in Seconds) :\n\n";
   lf << "Calculating Pressure Matrix : " << setprecision(2) << cpu_time[0].ClockTime() << endl;
   lf << "Solving Momentum Equations :  " << setprecision(2) << cpu_time[1].ClockTime();
   lf << " (" << cpu_time[1].ClockTime()/double(step) << " per time step)\n";
   lf << "Solving Pressure Equation :   " << setprecision(2) << cpu_time[2].ClockTime();
   lf << " (" << cpu_time[2].ClockTime()/double(step) << " per time step)\n";
   lf << "Updating Velocity :           " << setprecision(2) << cpu_time[3].ClockTime();
   lf << " (" << cpu_time[3].ClockTime()/double(step) << " per time step)\n";
   double tt = cpu_time[0].ClockTime() + cpu_time[1].ClockTime() + 
               cpu_time[2].ClockTime() + cpu_time[3].ClockTime();
   lf << "\nTOTAL TIME :   " << setprecision(2) << tt;

   delete Pp;
   if (vff) delete vff;
   if (pff) delete pff;
   if (sff) delete sff;
   if (wff) delete wff;
   if (vdf) delete vdf;
   lf.close();
   return 0;
}