weno/applications/euler/2d/VortexRotGFM/src/Problem.h

// -*- C++ -*-

// Copyright (C) 2003-2007 California Institute of Technology
// Ralf Deiterding, ralf@amroc.net

#ifndef AMROC_PROBLEM_H
#define AMROC_PROBLEM_H

#define DIM  2
  
#include "WENOProblem.h"

#define f_ipvel FORTRAN_NAME(ipvel, IPVEL)

extern "C" {
  void f_ipvel(); 
}

#define OWN_GFMAMRSOLVER
#include "WENOStdGFMProblem.h"
#include "F77Interfaces/F77GFMExactSolution.h"

class SolverSpecific : 
  public AMRGFMSolver<VectorType,FixupType,FlagType,DIM> {
  typedef VectorType::InternalDataType DataType;
  typedef AMRGFMSolver<VectorType,FixupType,FlagType,DIM> base;
  typedef GhostFluidMethod<VectorType,DIM>::boundary_conditions_type gfm_bc_type;
public:
  SolverSpecific(IntegratorSpecific& integ, 
                 base::initial_condition_type& init,
                 base::boundary_conditions_type& bc) : base(integ, init, bc) {
    SetLevelTransfer(new F77LevelTransfer<VectorType,DIM>(f_prolong, f_restrict));
#ifdef f_flgout
    SetFileOutput(new F77GFMFileOutput<VectorType,FixupType,FlagType,DIM>(*this,f_flgout)); 
#else   
    SetFileOutput(new GFMFileOutput<VectorType,FixupType,FlagType,DIM>(*this)); 
#endif
    SetFixup(new FixupSpecific());
    SetFlagging(new FlaggingSpecific(*this));   
    AddGFM(new GhostFluidMethod<VectorType,DIM>(
              new F77GFMBoundary<VectorType,DIM>(f_ibndrfl,f_itrans,f_ipvel),
              new F77GFMLevelSet<DataType,DIM>(f_lset),
              (gfm_bc_type*)(new F77BoundaryConditions<Vector<DataType,1>,DIM>(f_boundary))));
    SetExactSolution(new F77GFMExactSolution<VectorType,FixupType,FlagType,DIM>(*this,f_exact)); 
  }  
 
  ~SolverSpecific() {
    delete _GFM[0]->GetBoundaryConditionsP();
    DeleteGFM(_GFM[0]);
    delete _Flagging;
    delete _Fixup;
  }

  virtual void SetupData() {
    base::SetupData();
    _GFM[0]->BoundaryConditions_().SetupData(base::PGH(), base::NGhosts());
  }

  virtual void Initialize_(const double& dt_start) {   
    base::Initialize_(dt_start);
    CalculateCurrentMass(mass_old);
  }

  virtual double Tick(int VariableTimeStepping, const double dtv[], 
                      const double cflv[], int& Rejections) {

    double cfl = base::Tick(VariableTimeStepping,dtv,cflv,Rejections);

    VectorType mass_new;
    CalculateCurrentMass(mass_new);

    // Append to output file only on processor VizServer
    int me = MY_PROC; 
    if (me == VizServer) {
      std::ofstream outfile;
      std::ostream* out;
      std::string fname = "mass.dat";
      outfile.open(fname.c_str(), std::ios::out | std::ios::app);
      out = new std::ostream(outfile.rdbuf());
      *out << base::t[0];
      for (int n=0; n<base::NEquations(); n++)
        *out << " " << mass_old[n]-mass_new[n];
      *out << std::endl;
      outfile.close();
      delete out;      
    }

    return cfl;
  }

  void CalculateCurrentMass(VectorType& mass) {
    int Level=0;
    int Time = CurrentTime(base::GH(),Level); 
    int TStep = TimeStep(base::U(),Level);      
    forall (base::U(),Time,Level,c)   
      base::U()(Time+TStep,Level,c).equals(base::U()(Time,Level,c));
    end_forall

    if (base::NGFM()>0) {
      forall (base::U(),Time+TStep,Level,c)           
        for (int g=0; g<base::NGFM(); g++) {
          if (!base::GFM(g).IsUsed()) continue;
            BeginFastIndex2(phi, base::GFM(g).Phi()(Time,Level,c).bbox(), 
                            base::GFM(g).Phi()(Time,Level,c).data(), DataType);
            BeginFastIndex2(u, base::U()(Time+TStep,Level,c).bbox(), 
                            base::U()(Time+TStep,Level,c).data(), VectorType);
            for_2 (n, m, base::GFM(g).Phi()(Time,Level,c).bbox(), 
                   base::GFM(g).Phi()(Time,Level,c).bbox().stepsize())
              if (FastIndex2(phi,n,m)<-base::GFM(g).Boundary().Cutoff()) 
                FastIndex2(u,n,m) = static_cast<VectorType>(0.0);
            end_for
            EndFastIndex2(u);
            EndFastIndex2(phi);
          }
      end_forall
    }

    DCoords dx = base::GH().worldStep(base::U().GF_StepSize(Level));
    mass = Sum(base::U(),Time+TStep,Level);
    for (int d=0; d<base::Dim(); d++) mass *= dx(d);
  }

protected:
  VectorType mass_old;