sparse_grids_simulation.h

Go to the documentation of this file.

#ifndef _SPARSE_GRIDS_SIM_H_
#define _SPARSE_GRIDS_SIM_H_

#include <stdio.h>
#include <cstdlib>
#include <vector>
#include <string.h>
#include <string>
#include <utility>
#include <std_vec_ops.h>
#include <simulation.h>

#define _SPARSEGRIDS_SIM_INP_FNAME_ "sparse_grids.inp"

typedef std::vector<int> GridDimensions;
typedef std::vector<int> ProcDistribution;
typedef std::pair<double, GridDimensions> SGCTElem;

#define _coeff_ first
#define _dim_ second

class SparseGridsSimulation : public Simulation
{
  public:

    SparseGridsSimulation()
    {
      m_is_defined = false;

      m_nsims_sg = -1;
      m_ndims = -1;
      m_rank = -1;
      m_nproc = -1;

      m_sim_fg = NULL;
      m_sims_sg.clear();

      m_n_fg = -1;
      m_imin = -1;

      m_interp_order = 6;

      m_is_periodic.clear();
      m_combination.clear();
    }

    ~SparseGridsSimulation()
    {
      int err;
      /* clean up full grid object */
      err = CleanupBarebones(m_sim_fg);
      delete m_sim_fg;
      /* clean up and delete sparse grids objects */
      err = Cleanup((void*) m_sims_sg.data(), m_nsims_sg);
      if (err) {
        printf( "Error: CleanUp() returned with status %d on process %d.\n",
                err, m_rank );
      }
      m_sims_sg.clear();
      /* done */
      return;
    }

    int define(int, int);

    inline int ReadInputs()
    {
      int retval = ::ReadInputs(  (void*) m_sim_fg,
                                  1, 
                                  m_rank );
      return retval;
    }

    int Initialize();

    int InitialSolution();

    inline int InitializeBoundaries()
    {
      int retval = ::InitializeBoundaries(  (void*) m_sims_sg.data(),
                                            m_nsims_sg );
      if (m_nsims_sg > 0) {
        m_is_periodic.resize(m_ndims);
        for (int d=0; d<m_ndims; d++) {
          m_is_periodic[d] = ( m_sims_sg[0].solver.isPeriodic[d] == 1 ? true : false);
          m_sim_fg->solver.isPeriodic[d] = m_sims_sg[0].solver.isPeriodic[d];
        }
      }
      return retval;
    }

    inline int InitializeImmersedBoundaries()
    {
      int retval = ::InitializeImmersedBoundaries(  (void*) m_sims_sg.data(),
                                                    m_nsims_sg );
      return retval;
    }

    inline int InitializePhysics()
    {
      int retval = ::InitializePhysics( (void*) m_sims_sg.data(),
                                        m_nsims_sg );
      return retval;
    }

    inline int InitializePhysicsData()
    {
      for (int ns = 0; ns < m_sims_sg.size(); ns++) {
        int retval = ::InitializePhysicsData( (void*) &(m_sims_sg[ns]),
                                              -1, 
                                              m_nsims_sg,
                                              m_sim_fg->solver.dim_global);
        if (retval) {
          fprintf(stderr, "Error in SparseGridsSimulation::InitializePhysicsData()\n");
          fprintf(stderr, "  InitializePhysicsData returned with error code %d on rank %d.\n",
                  retval, m_sims_sg[ns].mpi.rank);
          return retval;
        }
      }
      return 0;
    }

    inline int InitializeSolvers()
    {
      int retval = ::InitializeSolvers( (void*) m_sims_sg.data(),
                                        m_nsims_sg );
      InitializeSolversBarebones(m_sim_fg);

      /* some modifications to output filename roots */
      for (int i=0; i<m_nsims_sg; i++) {
        strcpy(m_sims_sg[i].solver.op_fname_root, "op_sg");
        strcpy(m_sims_sg[i].solver.aux_op_fname_root, "ts0_sg");
      }
      strcpy(m_sim_fg->solver.op_fname_root, "op_fg");
      strcpy(m_sim_fg->solver.aux_op_fname_root, "ts0_fg");
      return retval;
    }

    int InitializationWrapup();

    int Solve();

    void WriteErrors(double, double);

    inline bool isDefined() const { return m_is_defined; }

#ifndef serial

    inline int mpiCommDup() 
    {
      if (!m_sim_fg) {
        fprintf(stderr, "Error in SparseGridsSimulation::mpiCommDup()\n");
        fprintf(stderr, "  m_sim_fg is not allocated on rank %d!\n", m_rank);
        return 1;
      }
      MPI_Comm_dup(MPI_COMM_WORLD, &(m_sim_fg->mpi.world));
      return 0;
    }
#endif

#ifdef with_petsc

    inline void usePetscTS(PetscBool a_flag)
    {
      m_sim_fg->solver.use_petscTS = a_flag;
    }

    inline int SolvePETSc()
    {
      int retval = ::SolvePETSc(  (void*) m_sims_sg.data(),
                                  m_nsims_sg,
                                  m_rank,
                                  m_nproc );
      return retval;
    }
#endif

  protected:

    bool  m_is_defined;   
    int   m_nsims_sg;     
    int   m_ndims;        
    int   m_rank,         
          m_nproc;        
    int   m_interp_order; 

    std::vector<bool> m_is_periodic; 
    int m_write_sg_solutions; 

    int m_print_sg_errors;

    SimulationObject* m_sim_fg;               
    std::vector<SimulationObject> m_sims_sg;  
    int m_n_fg; 
    int m_imin; 
    std::vector<SGCTElem> m_combination;  
    std::vector<ProcDistribution> m_iprocs; 
    int SanityChecks();

    int ComputeSGDimsAndCoeffs();

    int ComputeProcessorDistribution( ProcDistribution&,
                                      const GridDimensions&);

    void GetCTGridSizes(const int, std::vector<GridDimensions>&);

    int InitializeBarebones(SimulationObject*);

    int InitializeSolversBarebones(SimulationObject*);

    int CleanupBarebones(SimulationObject*);

    int SetSolverParameters(  SimulationObject&,
                              const GridDimensions&,
                              const ProcDistribution&,
                              const SimulationObject&,
                              const int,
                              const int );

    void OutputSolution(double);

    void CombinationTechnique(SimulationObject* const);

    void interpolate( SimulationObject* const, 
                      const SimulationObject* const);

    void interpolate( const GridDimensions&,
                      double** const, 
                      const SimulationObject* const);

    void interpolateGrid( SimulationObject* const, 
                          const SimulationObject* const);

    void coarsenGrid1D( const GridDimensions&,
                        const GridDimensions&,
                        const double* const, 
                        double* const,
                        int );

    void refineGrid1D( const GridDimensions&,
                       const GridDimensions&,
                       const double* const, 
                       double* const,
                       int );

    void fillGhostCells(  const GridDimensions&,
                          const int,
                          double* const,
                          const int );

    void CalculateError();

    void computeError( SimulationObject&, double* const);

    inline bool isPowerOfTwo(int x) 
    {
      if (x == 0)  return false;

      while (x > 1) {
        if (x%2 != 0) return false;
        x /= 2;
      }
      return true;
    }

    inline double cfunc(int a, int b)
    {
      return ( ( (double)factorial(a) ) / ( (double) (factorial(b)*factorial(a-b)) ) );
    }

    inline int factorial(int a)
    {
      int retval = 1.0;
      for (int i=1; i<=a; i++) retval *= i;
      return retval;
    }

    inline void allocateGridArrays( const GridDimensions& a_dim,
                                    double** const        a_x,
                                    const int             a_ngpt = 0)
    {
      GridDimensions dim_wghosts = a_dim;
      for (int i=0; i<dim_wghosts.size(); i++) {
        dim_wghosts[i] += (2*a_ngpt);
      }
      long size_x = StdVecOps::sum(dim_wghosts);
      (*a_x) = (double*) calloc (size_x, sizeof(double));
      for (int i=0; i<size_x; i++) (*a_x)[i] = 0.0;
      return;
    }

    inline void allocateDataArrays( const GridDimensions& a_dim,
                                    const int             a_nvars,
                                    double** const        a_u,
                                    const int             a_ngpt = 0)
    {
      GridDimensions dim_wghosts = a_dim;
      for (int i=0; i<dim_wghosts.size(); i++) {
        dim_wghosts[i] += (2*a_ngpt);
      }
      long size_u = a_nvars*StdVecOps::product(dim_wghosts);
      (*a_u) = (double*) calloc (size_u, sizeof(double));
      for (int i=0; i<size_u; i++) (*a_u)[i] = 0.0;
      return;
    }

  private:

};

#endif