#include "ionosphere.h"
#include "ionosphere_functions.h"
#include "DateTime.h"

/**
 * \file ionosphere.C
 * \brief The main program file.
 *
 * \mainpage 
 * \author Slava Merkin (vgm at bu.edu)
 * \since 2005
 * \date 24 August 2007
 *
 * The MIX code's primary function is to provide the inner boundary
 * condition for global MHD models, but it serves in fact a much more
 * diverse purpose. It can utilize existing empirical models of
 * ionospheric conductance to obtain the electrostatic potential
 * solution (solver mode), or it can receive conductances from an
 * external ionosphere-thermosphere model and then compute the
 * potential (solver/coupler mode). In the future, it is planned that
 * this coupler will mediate communication between other geospace
 * component models as well, i.e. Rice Convection Model (RCM) and LFM
 * (Lyon-Fedder-Mobarry). The major design consideration in developing
 * this code was to adopt the loose coupling concept,
 * i.e. plug-and-play should be made as simple as possible. Two
 * software packages are utilized in order to achieve this goal. The
 * first is <A
 * HREF="http://www.cs.umd.edu/projects/hpsl/chaos/ResearchAreas/ic/">InterComm</A>
 * developed by Alan Sussman's group at the University of Maryland at
 * College Park. InterComm provides means for data communication (no
 * file dumps!) between different executables which can be
 * parallelized and run on different systems. The other sofware tool
 * is <A HREF="http://www.llnl.gov/casc/Overture/">Overture</A>, an
 * object-oriented framework for solving partial differential
 * equations on overset grids developed by Bill Henshaw at LLNL. Both
 * packages make extensive use of <A
 * HREF="http://www.llnl.gov/casc/Overture/henshaw/documentation/App/manual/manual.html">A++/P++</A>,
 * a set of C++ classes allowing vector array operations in a
 * serial/parallel environment.
 *
 * The MIX code is not so much an Overture code as it is an A++
 * code. Overture's functionality is used in two places only: to
 * interpolate between grids and to obtain the elliptic solve for the
 * ionospheric potential. The caveat is that Overture's overset grid
 * technology, something that the software was designed for, is not
 * very useful when one needs to interpolate between entire grids (the
 * problem that we are dealing with). The problem that Overture is
 * solving is, in a way, much more sofisticated - to compute
 * interpolation boundaries. Fortunately, Overture does have simpler
 * interpolation functions, i.e. to interpolate from a grid onto a set
 * of unstructured points, and this is exactly what we are using. Any
 * interpolation sofware could be used for this purpose and we could
 * in principle use our own interpolation functions too. The second
 * place where we use Overture is where its use is more justified. The
 * application of Overture to obtain the ionospheric elliptic solution
 * makes this code very flexible and potentially gives nice
 * opportunities for interesting physical experiments: one can play
 * all kinds of games with different boundary conditions, with the
 * location of the low-latitude boundary, the grid resolution, etc. In
 * addition, Overture can interact with a number of sparse matrix
 * solvers (Yale, Harwell, SLAP, and PETSc) and a variety of matrix
 * inversion algorithms in combination with different preconditioners
 * can be used. An option unexplored so far, but a very attractive
 * one, is the use of multigrid solvers which are also implemented in
 * Overture.
 *
 * \image html structure.png
 * \image latex structure.eps "" width=\textwidth
 *
 * The picture above demonstrates the general structure of MIX and its
 * communication with other components. The box on the left represents
 * an unspecified global MHD model. The P++ interface in the MHD code
 * is not really a part of the MIX code (maybe it should be such) but
 * it is what actually talks to the MHD model. This interface
 * calculates the field-aligned current, the density and the sound
 * speed of plasma at the MHD inner boundary and sends these via
 * InterComm communication channel to the MIX code. After the
 * ionospheric potential solution is obtained, it is sent back to the
 * MHD interface which calculates plasma transverse velocity and
 * electric field and gives it back to the MHD code. The only
 * requirement to the MHD code so far is that its data be defined on a
 * logically spherical grid, but it is anticipated that an Overture
 * module will be developed to reduce this condition to an arbitrary
 * grid. The interface is a fully parallel P++ code, and thanks to the
 * InterComm beauty, the MIX code does not care about it. The
 * structure of the MIX code (box on the right) is similar. It has an
 * interface (actually, more than one - see the rest of the
 * documentation) which talks to an MHD or an ITM model (and maybe
 * something like RCM too), does all the data massaging, conversion of
 * physical variables, interpolation between different grids, etc. If
 * needed the ionospheric conductances may be calculated locally using
 * an empirical model (the conventional LFM
 * faxed-memo-kinda-AMIE-model or the Moen-Brekke, 1993 can be
 * used). Then all the data required are passed to the elliptic solver
 * which does the solve, hands the solution back to the interface, the
 * data required is extracted and is sent back to the components that
 * need it. Simple as that... Enjoy!
 *
 * <em>This code development was supported by the Center for
 * Integrated Space Weather Model (CISM) funded by the National
 * Science Foundation STC program under agreement ATM-012950.</em>
 */

/********************************************************************//**
 * \author Slava Merkin (vgm at bu.edu)
 * \since 2005 
 *
 * \param[in] argv[1] The name of the parameter file
 * 
 * This is the main program. There is nothing that should be brought
 * out from the source code to the manual; the comments in the code
 * are self-explanatory.
 ************************************************************************/
int main(int argc, char *argv[])  {
  DateTime currentUT;
  real tilt;
  clock_t startTime;

  Overture::start(argc,argv);  // initialize Overture

  // ===== Read MIX input parameters ========================================================
  char* paramFileName = argv[1];
  if ( !paramFileName ) {
    cerr << "Name of parameter file not entered. Usage: ionosphere <name of parameter file>" << endl;
    exit(1);
  } 
  else{
    cout << "Name of parameter file " << paramFileName << endl;
  }

  /* Read the parameter file and store parameter values in class */
  Params params(paramFileName);
  // MIX-ITM should start talking at
  DateTime itmCouplingStart = params.itmCouplingStart;
  // MIX-ITM should talk every "itmCouplingDelta" seconds;
  double itmCouplingDelta = params.itmCouplingDelta;
  // FIXME: Should spit out cerr warning message about LFM-MIX dump cadence versus itmCouplingDelta.
  
  // ===== Setup MHD & ITM interfaces =======================================================
  // Setup the interface to MHD & ITM codes
  MHD_Interface *mhd = NULL; 
  ITM_Interface *itm = NULL;
  IM_Interface  *im  = NULL;

  int rank = 0;
  int ic_err;
  IC_XJD* xjd;

  if (params.use_intercomm){
    // Initialize InterComm    
    xjd = IC_Initialize(params.localName,rank,params.xjdName,&ic_err);
    std::cout << "Initialized with " << params.localName << " " << params.xjdName << "\n";
    mhd = new MHD_IC_Interface(&params, xjd);
    im  = new IM_IC_Interface(xjd);
    itm = new ITM_IC_Interface(&params, xjd);
    
    // We commit IC regions here because this function commits all
    // regions created previously.  In our case, the MHD,
    // Ionosphere-Theremosphere-Mesosphere and Inner-Magnetosphere
    // interfaces create IC_regions that should be commited.
    IC_Commit_region(xjd,&ic_err);

    // Inner Magnetosphere must initialize after Commit_region & before main loop
    im->Initialize();
  }
  else{
    // Initialize file exchanges
    mhd = new MHD_FE_Interface(&params, xjd);
    itm = new ITM_FE_Interface(&params,xjd);
    cerr << "***** ERROR ***** RCM coupling requires the InterComm interface, but you have selected file exchange!" << endl;
  }

  // ===== Initialize MIX Grids =============================================================   
  /* Call the Solver constructor The parameters are defined in const.h. See
   * comments in that file for requirements */
  Solver ionSolver(params.mixGrid.getNLAT(),
		   params.mixGrid.getLATITUDE(),
		   params.mixGrid.getSOLVER_LON_GRID_STEP() );

  // Are we coupling to an ITM model?
  CompositeGrid smGlobe;
  realCompositeGridFunction smGlobeFunctions;
  if (itmCouplingStart.getYear() < 9998){
    /* 
     *  These are created on 2007-15-08 in order to test the following:
     *  When stuff is interpolated from the solver grid onto the ITM
     *  grid the points at mid- and low latitudes (those in between the
     *  southern and northern solver grids) are extrapolated from the
     *  boundary. My idea was to create a grid (theta, phi) covering the
     *  entire globe in SM coordinates and set these points to zero (see
     *  function Solver::mapToSMGlobe called below). After testing this
     *  possibility, I am now not sure whether to go this way, because
     *  it may create sharp gradients at the boundary. Need to check
     *  with Wenbin and Mike. Meanwhile, I'm going to do it this way
     *  because it really simplifies the ITM interface (interpolation is
     *  now done in spherical angular coordinates and there is no
     *  division into hemisphere. In fact, I have removed ITM_Interface
     *  interpolation functions based on the cartesian grids.
     */
    
    smGlobe = ionSolver.generateSMGlobeGrid();
    //  realCompositeGridFunction psi(smGlobe), averageEnergy(smGlobe),numberFlux(smGlobe);
    
    //  smGlobeFunctions are stored as different components:
    //  0: potential;  1: average energy;  2: number flux (see functions.h)
    smGlobeFunctions = realCompositeGridFunction(smGlobe,all,all,all,3);
  }

  // Construct Solver grid
  CompositeGrid sgrid = ionSolver.generateCartesian2DGrid();
  
  // Get the second shell MHD ionospheric grid from the MHD Interface
  CompositeGrid cg = mhd->getSecondShellIonCartesianGrid();
  realCompositeGridFunction secondShellMHDFunctions(cg,all,all,all,8);

  // Get the first shell MHD ionospheric grid from the MHD interface
  CompositeGrid cg1 = mhd->getFirstShellIonCartesianGrid();
  realCompositeGridFunction firstShellMHDFunctions(cg1,all,all,all,2);

  realCompositeGridFunction solverFunctions(sgrid,all,all,all,18);
  solverFunctions.setName("Solver grid functions");
  solverFunctions.setName("Potential North [V]",SOLVER_POTENTIAL_NORTH);         // MIX: phi
  solverFunctions.setName("Potential South [V]",SOLVER_POTENTIAL_SOUTH);         //      phi
  solverFunctions.setName("FAC North [A/m^2]",SOLVER_FAC_NORTH);                 // ITM: J_w
  solverFunctions.setName("FAC South [A/m^2]",SOLVER_FAC_SOUTH);                 //      J_w
  solverFunctions.setName("Sound speed North [cm/s]",SOLVER_SOUND_SPEED_NORTH);  // MHD: c_s
  solverFunctions.setName("Sound speed South [cm/s]",SOLVER_SOUND_SPEED_SOUTH);  //      c_s
  solverFunctions.setName("Density North [g/cm^3]",SOLVER_DENSITY_NORTH);        // MHD: rho
  solverFunctions.setName("Density South [g/cm^3]",SOLVER_DENSITY_SOUTH);        //      rho
  solverFunctions.setName("Pedersen conductance North [S]",SOLVER_SIGMAP_NORTH); // ITM: sigma_p
  solverFunctions.setName("Pedersen conductance South [S]",SOLVER_SIGMAP_SOUTH); //      sigma_p
  solverFunctions.setName("Hall conductance North [S]",SOLVER_SIGMAH_NORTH);     // ITM: sigma_h
  solverFunctions.setName("Hall conductance South [S]",SOLVER_SIGMAH_SOUTH);     //      sigma_h
  solverFunctions.setName("Average energy North [keV]",SOLVER_AVG_ENG_NORTH);    // MIX: E
  solverFunctions.setName("Average energy South [keV]",SOLVER_AVG_ENG_SOUTH);    //      E
  solverFunctions.setName("Number flux North [1/cm^2 s]",SOLVER_NUM_FLUX_NORTH); // MIX: F
  solverFunctions.setName("Number flux South [1/cm^2 s]",SOLVER_NUM_FLUX_SOUTH); //      F
  solverFunctions.setName("Neutral wind speed North [m/s]",SOLVER_WIND_NORTH);   //
  solverFunctions.setName("Neutral wind speed South [m/s]",SOLVER_WIND_SOUTH);   //
  
  // Constructor for the conductance model
  // If an ITM is used to calculate the conductances, then
  // all the constants are just ignored: all we need is the solver info
  Conductance sigma(&ionSolver,&params);

  // ===== Setup IO/Overture visualization ==================================================
#ifdef DX_OUT
  IO DXHDF(&params);
#else
  Ogshow show(params.outputFilePrefix,".",1);
  show.setFlushFrequency(params.flushFrequency);
#endif

  // ===== MHD Interface Helpers ============================================================
  unsigned long int mhdTimeStep;
  unsigned long int mhdExchangeCount=1; // LFM starts counting loops from 1, so MIX should too.
  bool isLastStep;

  // ===== ITM Interface helpers ============================================================
  DateTime itmCouplingNext = itmCouplingStart;
  unsigned long int itmExchangeCount = 0;  // How many MIX-ITM exchanges?

  // ===== IM Interface helpers =============================================================
  int exchangeWithIm = 0;
  unsigned long int imExchangeCount = 0;

  // Timekeeping things:
  double t_mhd2MIX, t_FandE, t_conduct, t_solve, t_MIX2mhd, t_MIX2itm, t_MIX2im, t_IO;
  
  // ===== MAIN LOOP ========================================================================
  // The loop will continue until KILL signal received from the MHD code
  while ( receiveMHDscalars(mhd, currentUT, tilt, mhdTimeStep, isLastStep, exchangeWithIm, imExchangeCount) )
    {
      std::cout << "LFM-MIX exchange # " << mhdExchangeCount << ";  "
		<< "MHD UT: " <<  currentUT.getDateTimeString() << endl;

      t_mhd2MIX  = ionosphere_MHD_to_MIX(mhd,secondShellMHDFunctions, solverFunctions, ionSolver); 

      t_FandE    = ionosphere_get_Flux_and_Energy(exchangeWithIm, im, sigma, solverFunctions, tilt);

      t_conduct  = ionosphere_get_Conductances(params, mhdExchangeCount, itmExchangeCount, currentUT, itmCouplingNext,
      					       sigma, solverFunctions, DXHDF, tilt, itm);      

      t_solve    = ionosphere_compute_Potential(ionSolver, solverFunctions);

      t_MIX2mhd  = ionosphere_MIX_to_MHD(mhd,ionSolver, solverFunctions, firstShellMHDFunctions,secondShellMHDFunctions);

      t_MIX2itm  = ionosphere_MIX_to_ITM(params, isLastStep, itmExchangeCount, currentUT, itmCouplingNext, itm,
					 ionSolver, solverFunctions, smGlobeFunctions, DXHDF, tilt);

      t_MIX2im   = ionosphere_MIX_to_IM(exchangeWithIm, im, solverFunctions);
      
      startTime = clock();
#ifdef DX_OUT
      if ( ( (mhdExchangeCount-1) % params.flushFrequency ) == 0 ){	
	DXHDF.Write(tilt, currentUT, mhdTimeStep, solverFunctions);
      }
#else
      show.startFrame();
      char label[STRLEN];
      
      int nsym=sprintf(label,"%s%08d\n","MHD time step: ",mhdTimeStep);
      show.saveComment(0,label);
      nsym=sprintf(label,"%04d//%02d//%02d %02d:%02d:%02f\n",
		   currentUT.getYear(), currentUT.getMonth(),  currentUT.getDay(),
		   currentUT.getHours(),currentUT.getMinutes(),currentUT.getSeconds());
      show.saveComment(1,label);

      show.saveSolution(solverFunctions);
      show.endFrame();
#endif // DX_OUT
      t_IO = (float) (clock()-startTime)/CLOCKS_PER_SEC;

      // Spit out timing information for this step:
      std::cout << "\t================================================================================\n"
		<< "\tdt = " << (t_mhd2MIX + t_FandE + t_conduct + t_solve + t_MIX2mhd + t_MIX2itm + t_MIX2im + t_IO) << "; "
		<< " (MHD->MIX interp: " << t_mhd2MIX
		<< ", Flux & Energy: " << t_FandE
		<< ", Conductances: " << t_conduct
		<< ", MIX solve: " << t_solve
		<< ", MIX->MHD: " << t_MIX2mhd 
		<< ", MIX->ITM: " << t_MIX2itm
		<< ", MIX->IM: " << t_MIX2im
		<< ", IO: " << t_IO
		<< ")\n"
		<< "\t================================================================================" << endl;
      
      mhdExchangeCount+=1;

    }
  // ===== MAIN LOOP ENDS HERE ==============================================================
  // ========================================================================================

  delete itm;
  itm = NULL;

  delete mhd;
  mhd = NULL;

  delete im;
  im = NULL;

  if (params.use_intercomm){
    IC_Finalize(xjd,&ic_err);
  }

  Overture::finish();
  return 0;
}