gsXBraid.h

Go to the documentation of this file.

#pragma once

#include <gismo.h>

#ifndef GISMO_WITH_MPI
#define braid_SEQUENTIAL 1
#endif

#include <XBraid/braid/braid.hpp>

namespace gismo {

  class gsXBraidAccessStatus;
  class gsXBraidSyncStatus;
  class gsXBraidStepStatus;
  class gsXBraidCoarsenRefStatus;
  class gsXBraidBufferStatus;
  class gsXBraidObjectiveStatus;
  
  template<typename T>
  class gsXBraid : public BraidApp
  {
  public:
    gsXBraid(const gsMpiComm& comm,
             const braid_Real tstart,
             const braid_Real tstop,
             braid_Int        ntime);
       
    virtual ~gsXBraid();

    virtual braid_Int Free(braid_Vector) { return braid_Int(0); }

    virtual braid_Int Residual(braid_Vector, braid_Vector, BraidStepStatus&)
    { GISMO_NO_IMPLEMENTATION }
    
    void solve() { core.Drive(); }
    
  public:
    void SetMaxLevels(braid_Int max_levels) { core.SetMaxLevels(max_levels); }
    
    void SetIncrMaxLevels() { core.SetIncrMaxLevels(); }

    void SetSkip(braid_Int skip) { core.SetSkip(skip); }

    void SetMinCoarse(braid_Int min_coarse) { core.SetMinCoarse(min_coarse); }

    void SetNRelax(braid_Int level, braid_Int nrelax) { core.SetNRelax(level, nrelax); }

    void SetNRelax(braid_Int nrelax) { core.SetNRelax(-1, nrelax); }
        
    void SetAbsTol(braid_Real tol) { core.SetAbsTol(tol); }

    void SetRelTol(braid_Real tol) { core.SetRelTol(tol); }

    void SetTemporalNorm(braid_Int tnorm) { core.SetTemporalNorm(tnorm); }

    void SetCFactor(braid_Int level, braid_Int cfactor) { core.SetCFactor(level, cfactor); }

    void SetCFactor(braid_Int cfactor) { core.SetCFactor(-1, cfactor); }

    void SetPeriodic(braid_Int periodic) { core.SetPeriodic(periodic); }

    void SetMaxIter(braid_Int max_iter) { core.SetMaxIter(max_iter); }

    void SetPrintLevel(braid_Int print_level) { core.SetPrintLevel(print_level); }

    void SetPrintFile(const char *printfile_name) { core.SetPrintFile(printfile_name); }
    
    void SetSeqSoln(braid_Int use_seq_soln) { core.SetSeqSoln(use_seq_soln); }

    void SetAccessLevel(braid_Int access_level) { core.SetAccessLevel(access_level); }

    void SetFMG() { core.SetFMG(); }

    void SetNFMG(braid_Int k) { core.SetNFMG(k); }

    void SetNFMGVcyc(braid_Int nfmg_Vcyc) { core.SetNFMGVcyc(nfmg_Vcyc); }

    void SetStorage(braid_Int storage) { core.SetStorage(storage); }

    void SetRefine(braid_Int refine) {core.SetRefine(refine);}

    void SetMaxRefinements(braid_Int max_refinements) {core.SetMaxRefinements(max_refinements);}

    void SetRichardsonEstimation(braid_Int est_error, braid_Int richardson, braid_Int local_order)
    { core.SetRichardsonEstimation(est_error, richardson, local_order); }

  public:
    void SetResidual() { core.SetResidual(); }
    
    void SetSpatialCoarsenAndRefine() { core.SetSpatialCoarsenAndRefine(); }
    
    void SetSync() { core.SetSync(); }

    void SetDefaultPrintFile() { core.SetDefaultPrintFile(); }

    void SetFileIOLevel(braid_Int io_level) { core.SetFileIOLevel(io_level); }

    void SetCRelaxWt(braid_Int level, braid_Real  Cwt) { core.SetCRelaxWt(level, Cwt); }

   void SetTPointsCutoff(braid_Int tpoints_cutoff) { core.SetTPointsCutoff(tpoints_cutoff); }

   void SetFullRNormRes(braid_PtFcnResidual residual) { core.SetFullRNormRes(residual); }

   void SetTimeGrid(braid_PtFcnTimeGrid tgrid) { core.SetTimeGrid(tgrid); }
    
  public:
    void GetNumIter(braid_Int *niter_ptr) { core.GetNumIter(niter_ptr); }

    void GetRNorms(braid_Int *nrequest_ptr, braid_Real *rnorms) { core.GetRNorms(nrequest_ptr, rnorms); }

    void GetNLevels(braid_Int *nlevels_ptr) { core.GetNLevels(nlevels_ptr); }

    void GetMyID(braid_Int *myid_ptr) { core.GetMyID(myid_ptr); }
    
    braid_Int iterations() {
      braid_Int niter;
      GetNumIter(&niter);
      return niter;
    }

    braid_Real norm(braid_Int nrequest) {
      braid_Real rnorm;
      GetRNorms(&nrequest, &rnorm);
      return rnorm;
    }
    
    braid_Int levels() {
      braid_Int nlevels;
      GetNLevels(&nlevels);
      return nlevels;
    }

    braid_Int id() {
      braid_Int myid;
      GetMyID(&myid);
      return myid;
    }
    
  protected:
    BraidCore core;    
  };

  
  template<typename T>
  class gsXBraid< gsMatrix<T> > : public gsXBraid<T>
  {
  public:
    gsXBraid(const gsMpiComm& comm,
             const braid_Real tstart,
             const braid_Real tstop,
             braid_Int        ntime);
    
    virtual ~gsXBraid();
    
    virtual braid_Int Clone(braid_Vector  u,
                            braid_Vector *v_ptr)
    {
      gsMatrix<T>* u_ptr = (gsMatrix<T>*) u;
      gsMatrix<T>* v     = new gsMatrix<T>();
      *v = *u_ptr;
      *v_ptr = (braid_Vector) v;
      return braid_Int(0);
    }
    
    virtual braid_Int Free(braid_Vector u)
    {
      gsMatrix<T>* u_ptr = (gsMatrix<T>*) u;
      delete u_ptr;
      return braid_Int(0);
    }
    
    virtual braid_Int Sum(braid_Real   alpha,
                          braid_Vector x,
                          braid_Real   beta,
                          braid_Vector y)
    {
      gsMatrix<T>* x_ptr = (gsMatrix<T>*) x;
      gsMatrix<T>* y_ptr = (gsMatrix<T>*) y;
      *y_ptr = (T)alpha * (*x_ptr) + (T)beta * (*y_ptr);
      return braid_Int(0);
    }
    
    virtual braid_Int SpatialNorm(braid_Vector  u,
                                  braid_Real   *norm_ptr)
    {
      gsMatrix<T> *u_ptr = (gsMatrix<T>*) u;    
      *norm_ptr = u_ptr->norm();
      return braid_Int(0);
    }
    
    virtual braid_Int BufPack(braid_Vector       u,
                              void              *buffer,
                              BraidBufferStatus &status)
    {
      gsMatrix<T> *u_ptr = (gsMatrix<T>*) u;
      T* buffer_ptr      = (T*) buffer;
      T* data_ptr        = u_ptr->data();
      index_t size       = u_ptr->rows()*u_ptr->cols();
      
      buffer_ptr[0] = u_ptr->rows();
      buffer_ptr[1] = u_ptr->cols();
      for (index_t idx = 0; idx < size; ++idx)
        buffer_ptr[idx+2] = data_ptr[idx];
      
      status.SetSize(sizeof(T)*(size+2));
      return braid_Int(0);
    }

    virtual braid_Int BufUnpack(void              *buffer,
                                braid_Vector      *u_ptr,
                                BraidBufferStatus &status)
    {
      T* buffer_ptr  = (T*) buffer;
      index_t rows   = buffer_ptr[0];
      index_t cols   = buffer_ptr[1];
      gsMatrix<T>* u = new gsMatrix<T>(rows,cols);
      T* data_ptr       = u->data();

    for (index_t idx = 0; idx < rows*cols; ++idx)
      data_ptr[idx] = buffer_ptr[idx+2];
    
    *u_ptr = (braid_Vector) u;
    return braid_Int(0);
    }
  };

  template<typename T>
  class gsXBraid< gsVector<T> > : public gsXBraid<T>
  {
  public:
    gsXBraid(const gsMpiComm& comm,
             const braid_Real tstart,
             const braid_Real tstop,
             braid_Int        ntime);
    
    virtual ~gsXBraid();
    
    virtual braid_Int Clone(braid_Vector  u,
                            braid_Vector *v_ptr)
    {
      gsVector<T>* u_ptr = (gsVector<T>*) u;
      gsVector<T>* v     = new gsVector<T>();
      *v = *u_ptr;
      *v_ptr = (braid_Vector) v;
      return braid_Int(0);
    }
    
    virtual braid_Int Free(braid_Vector u)
    {
      gsVector<T>* u_ptr = (gsVector<T>*) u;
      delete u_ptr;
      return braid_Int(0);
    }
    
    virtual braid_Int Sum(braid_Real   alpha,
                          braid_Vector x,
                          braid_Real   beta,
                          braid_Vector y)
    {
      gsVector<T>* x_ptr = (gsVector<T>*) x;
      gsVector<T>* y_ptr = (gsVector<T>*) y;
      *y_ptr = (T)alpha * (*x_ptr) + (T)beta * (*y_ptr);
      return braid_Int(0);
    }
    
    virtual braid_Int SpatialNorm(braid_Vector  u,
                                  braid_Real   *norm_ptr)
    {
      gsVector<T> *u_ptr = (gsVector<T>*) u;    
      *norm_ptr = u_ptr->norm();
      return braid_Int(0);
    }
    
    virtual braid_Int BufPack(braid_Vector       u,
                              void              *buffer,
                              BraidBufferStatus &status)
    {
      gsVector<T> *u_ptr = (gsVector<T>*) u;
      T* buffer_ptr      = (T*) buffer;
      T* data_ptr        = u_ptr->data();
      index_t size       = u_ptr->size();
      
      buffer_ptr[0] = u_ptr->size();
      for (index_t idx = 0; idx < size; ++idx)
        buffer_ptr[idx+1] = data_ptr[idx];
      
      status.SetSize(sizeof(T)*(size+1));
      return braid_Int(0);
    }

    virtual braid_Int BufUnpack(void              *buffer,
                                braid_Vector      *u_ptr,
                                BraidBufferStatus &status)
    {
      T* buffer_ptr  = (T*) buffer;
      index_t size   = buffer_ptr[0];
      gsVector<T>* u = new gsVector<T>(size);
      T* data_ptr    = u->data();

    for (index_t idx = 0; idx < size; ++idx)
      data_ptr[idx] = buffer_ptr[idx+1];
    
    *u_ptr = (braid_Vector) u;
    return braid_Int(0);
    }
  };
  
  class gsXBraidAccessStatus : public BraidAccessStatus
  {
  public:    
    braid_Int iterations() {
      braid_Int iter;
      GetIter(&iter);
      return iter;
    }

    braid_Int level() {
      braid_Int level;
      GetLevel(&level);
      return level;
    }

    braid_Int levels() {
      braid_Int nlevels;
      GetNLevels(&nlevels);
      return nlevels;
    }

    braid_Int refines() {
      braid_Int nref;
      GetNRefine(&nref);
      return nref;
    }

    braid_Real time() {
      braid_Real t;
      GetT(&t);
      return t;
    }
    
    braid_Int times() {
      braid_Int ntpoints;
      GetNTPoints(&ntpoints);
      return ntpoints;
    }
    
    bool done() {
      braid_Int status;
      GetDone(&status);
      return bool(status);
    }

    braid_Int callingFunction() {
      braid_Int callingfcn;
      GetCallingFunction(&callingfcn);
      return callingfcn;
    }
    
    braid_Int timeIndex() {
      braid_Int tindex;
      GetTIndex(&tindex);
      return tindex;
    }

    braid_Int test() {
      braid_Int wtest;
      GetWrapperTest(&wtest);
      return wtest;
    }

    braid_Real norm() {
      braid_Real rnorm;
      GetResidual(&rnorm);
      return rnorm;
    }

    braid_Real error() {
      braid_Real errorest;
      GetSingleErrorEstAccess(&errorest);
      return errorest;
    }
  };

  class gsXBraidSyncStatus : public BraidSyncStatus
  {
  public:
    braid_Int iterations() {
      braid_Int iter;
      GetIter(&iter);
      return iter;
    }

    braid_Int level() {
      braid_Int level;
      GetLevel(&level);
      return level;
    }

    braid_Int levels() {
      braid_Int nlevels;
      GetNLevels(&nlevels);
      return nlevels;
    }

    braid_Int refines() {
      braid_Int nref;
      GetNRefine(&nref);
      return nref;
    }
    
    braid_Int times() {
      braid_Int ntpoints;
      GetNTPoints(&ntpoints);
      return ntpoints;
    }
    
    bool done() {
      braid_Int status;
      GetDone(&status);
      return bool(status);
    }

    braid_Int callingFunction() {
      braid_Int callingfcn;
      GetCallingFunction(&callingfcn);
      return callingfcn;
    }

    braid_Real errors() {
      braid_Real errorest;
      GetAllErrorEst(&errorest);
      return errorest;
    }

    braid_Int nerrors() {
      braid_Int numerrorest;
      GetNumErrorEst(&numerrorest);
      return numerrorest;
    }
  };

  class gsXBraidStepStatus : public BraidStepStatus
  {
  public:
    braid_Int iterations() {
      braid_Int iter;
      GetIter(&iter);
      return iter;
    }

    braid_Int level() {
      braid_Int level;
      GetLevel(&level);
      return level;
    }

    braid_Int levels() {
      braid_Int nlevels;
      GetNLevels(&nlevels);
      return nlevels;
    }

    braid_Int refines() {
      braid_Int nref;
      GetNRefine(&nref);
      return nref;
    }

    braid_Real time() {
      braid_Real t;
      GetT(&t);
      return t;
    }
    
    braid_Int times() {
      braid_Int ntpoints;
      GetNTPoints(&ntpoints);
      return ntpoints;
    }

    braid_Real timeStop() {
      braid_Real t;
      GetTstop(&t);
      return t;
    }

    std::pair<braid_Real, braid_Real> timeInterval() {
      std::pair<braid_Real, braid_Real> t;
      GetTstartTstop(&t.first, &t.second);
      return t;
    }
    
    braid_Int timeIndex() {
      braid_Int tindex;
      GetTIndex(&tindex);
      return tindex;
    }

    braid_Real tol() {
      braid_Real t;
      GetTol(&t);
      return t;
    }

    braid_Real tolFine() {
      braid_Real t;
      GetOldFineTolx(&t);
      return t;
    }

    braid_Real error() {
      braid_Real errorest;
      GetSingleErrorEstStep(&errorest);
      return errorest;
    }

    braid_Real accuracy(braid_Real loose_tol, braid_Real tight_tol) {
      braid_Real tol;
      GetSpatialAccuracy(loose_tol, tight_tol, &tol);
      return tol;
    }
  };

  class gsXBraidCoarsenRefStatus : public BraidCoarsenRefStatus
  {
  public:
    braid_Int iterations() {
      braid_Int iter;
      GetIter(&iter);
      return iter;
    }

    braid_Int level() {
      braid_Int level;
      GetLevel(&level);
      return level;
    }

    braid_Int levels() {
      braid_Int nlevels;
      GetNLevels(&nlevels);
      return nlevels;
    }

    braid_Int refines() {
      braid_Int nref;
      GetNRefine(&nref);
      return nref;
    }

    braid_Real time() {
      braid_Real t;
      GetT(&t);
      return t;
    }
    
    braid_Int times() {
      braid_Int ntpoints;
      GetNTPoints(&ntpoints);
      return ntpoints;
    }

    braid_Int timeIndex() {
      braid_Int tindex;
      GetTIndex(&tindex);
      return tindex;
    }

    braid_Real ftimeStop() {
      braid_Real t;
      GetFTstop(&t);
      return t;
    }

    braid_Real ftimeStart() {
      braid_Real t;
      GetFTprior(&t);
      return t;
    }
    
    braid_Real ctimeStop() {
      braid_Real t;
      GetCTstop(&t);
      return t;
    }

    braid_Real ctimeStart() {
      braid_Real t;
      GetCTprior(&t);
      return t;
    }
  };

  class gsXBraidBufferStatus : public BraidBufferStatus
  {
  public:
    braid_Int type() {
      braid_Int msg;
      GetMessageType(&msg);
      return msg;
    }
  };

  class gsXBraidObjectiveStatus : public BraidObjectiveStatus
  {
  public:
    braid_Int iterations() {
      braid_Int iter;
      GetIter(&iter);
      return iter;
    }

    braid_Int level() {
      braid_Int level;
      GetLevel(&level);
      return level;
    }

    braid_Int levels() {
      braid_Int nlevels;
      GetNLevels(&nlevels);
      return nlevels;
    }

    braid_Int refines() {
      braid_Int nref;
      GetNRefine(&nref);
      return nref;
    }
    
    braid_Real time() {
      braid_Real t;
      GetT(&t);
      return t;
    }
    
    braid_Int times() {
      braid_Int ntpoints;
      GetNTPoints(&ntpoints);
      return ntpoints;
    }
    
    braid_Int timeIndex() {
      braid_Int tindex;
      GetTIndex(&tindex);
      return tindex;
    }
  };
  
}// namespace gismo

#ifndef GISMO_BUILD_LIB
#include GISMO_HPP_HEADER(gsXBraid.hpp)
#endif