clawpack/applications/scaling/2d/Ramp/src/Problem.h

// -*- C++ -*-

// Copyright (C) 2002 Ralf Deiterding
// Brandenburgische Universitaet Cottbus

#ifndef AMROC_PROBLEM_H
#define AMROC_PROBLEM_H

#include "euler2.h"

#undef  NAUX
#define NAUX      4

#include "ClpProblem.h"

#define OWN_FLAGGING
#define OWN_AMRSOLVER
#include "ClpStdProblem.h"

class ScaledGradientSpecific : public ScaledGradient<VectorType,FlagType,DIM> {
  typedef ScaledGradient<VectorType,FlagType,DIM> base;
  typedef VectorType::InternalDataType DataType;
public:
  typedef base::grid_fct_type grid_fct_type;
  typedef base::flag_fct_type flag_fct_type;

  virtual bool FlagByScaledGradient(grid_fct_type& work, flag_fct_type& flags,
                                    const int& Time, const int& Level, 
                                    const DataType& TolVal, const FlagType& FlagValue) { 
    if (TolVal <= 0.0) return false;
    int TStep = TimeStep(work,Level);
    forall(work,Time,Level,c)
      work(Time+TStep,Level,c) = 0.0;
    end_forall 
    int East[2]  = { 1, 0 }; int West[2]  = { -1, 0 };
    int North[2] = { 0, 1 }; int South[2] = { 0, -1 };
    base::Difference(work, Time, Level, East, West);
    base::Difference(work, Time, Level, North, South);      
    base::FlagByValue(work, flags, Time, Level, TolVal, FlagValue);
    return true;
  }
};

class FlaggingSpecific : 
  public AMRFlagging<VectorType,FixupType,FlagType,DIM> {
  typedef AMRFlagging<VectorType,FixupType,FlagType,DIM> base;
  typedef base::solver_type solver_type;
public:
  FlaggingSpecific(solver_type& solver) : base(solver) {
      base::AddCriterion(new ScaledGradientSpecific());
      base::AddCriterion(new AbsoluteError<VectorType,FixupType,FlagType,DIM>(solver));
    }
  ~FlaggingSpecific() { DeleteAllCriterions(); }
};

class SolverSpecific : 
  public AMRSolver<VectorType,FixupType,FlagType,DIM> {
  typedef AMRSolver<VectorType,FixupType,FlagType,DIM> base;
public:
  SolverSpecific(IntegratorSpecific& integ, 
                 base::initial_condition_type& init,
                 base::boundary_conditions_type& bc) :
    AMRSolver<VectorType,FixupType,FlagType,DIM>(integ, init, bc) {
    SetLevelTransfer(new F77LevelTransfer<VectorType,DIM>(f_prolong, f_restrict));
#ifdef f_flgout
    SetFileOutput(new F77FileOutput<VectorType,DIM>(f_flgout)); 
#else   
    SetFileOutput(new FileOutput<VectorType,DIM>()); 
#endif
    SetFixup(new FixupSpecific(integ));
    SetFlagging(new FlaggingSpecific(*this)); 
    last[0]=0.; last[1]=0.; last[2]=0.;
    TraceDumpEvery = 1;
  }  

  ~SolverSpecific() {
    delete _LevelTransfer;
    delete _Flagging;
    delete _Fixup;
    delete _FileOutput;
  }

  virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
    base::register_at(Ctrl,prefix);
    RegisterAt(base::LocCtrl,"TraceDumpEvery",TraceDumpEvery);
  } 
  virtual void register_at(ControlDevice& Ctrl) {
    base::register_at(Ctrl);
  }

  virtual void Initialize(double& t, double& dt) {
    base::Initialize(t,dt);
    TraceDump(0);
  }

  virtual void Advance(double& t, double& dt) {
    base::Advance(t,dt);
    int TimeZero  = CurrentTime(base::GH(),0);
    TimeZero /= RefinedBy(base::GH(),MaxLevel(base::GH()));
    if (TimeZero%TraceDumpEvery==0 && TraceDumpEvery>0) 
      TraceDump(TimeZero/TraceDumpEvery);
  }

  void TraceDump(int time) {
    int me = MY_PROC; 
#ifdef TIMING_AMR
    double stats[3], current[3];
    timing.collect_timing(comm_service::comm(),Timing::INTEGRATION_MAIN,current);
    timing.collect_timing(comm_service::comm(),Timing::INTEGRATION_ESTIMATE,stats);
    if (me == VizServer) {
      current[0]+=stats[0]; current[1]+=stats[1]; current[2]+=stats[2];
    }
    timing.collect_timing(comm_service::comm(),Timing::INTEGRATION_SHADOW,stats);
    if (me == VizServer) {
      current[0]+=stats[0]; current[1]+=stats[1]; current[2]+=stats[2];
    }
    timing.collect_timing(comm_service::comm(),Timing::SOURCE_INTEGRATION,stats);
    if (me == VizServer) {
      current[0]+=stats[0]; current[1]+=stats[1]; current[2]+=stats[2];
    }
    if (me == VizServer) {
      stats[0] = current[0]-last[0]; stats[1] = current[1]-last[1]; stats[2] = current[2]-last[2]; 
      std::ofstream ofst("balance.txt", std::ios::out | std::ios::app);
      char str[128];
      std::sprintf(str, "\t%3.9f    %2.4f    %2.4f    %2.4f", base::t[0], 
                   (stats[2]>0. ? stats[0]/stats[2] : 1.),
                   (stats[2]>0. ? stats[1]/stats[2] : 1.),
                   (stats[2]>0. ? (stats[0]-stats[1])/stats[2] : 1.));
      ofst << time << str << std::endl;
      ofst.close();
    }
#endif

    if (me == VizServer) {
      int id=0;
      std::ofstream ofs("trace.txt", std::ios::out | std::ios::app);
      for (register int l=0; l<=FineLevel(base::GH()); l++)
        for (GridBox *gb=base::GH().ggbl(l)->first();gb;gb=base::GH().ggbl(l)->next())
          dumpBox(gb->gbBBox(), l, gb->gbOwner(), time, id++, ofs);
      ofs.close();
    }
  }

  void dumpBox(const BBox &bb, int l, int proc, int time, int id, std::ofstream &ofs) {
    ofs << time << "\t" << proc << "\t" << l << "\t" << id <<"\t";
    ofs << bb.extents(0) <<"\t";
    if (bb.rank>1) ofs << bb.extents(1) <<"\t";
    else ofs << "-1\t";
    if (bb.rank>2) ofs << bb.extents(2) <<"\t";
    else ofs << "-1\t";
    ofs << bb.lower(0) <<"\t";
    if (bb.rank>1) ofs << bb.lower(1) <<"\t";
    else ofs << "-1\t";
    if (bb.rank>2) ofs << bb.lower(2) <<"\t";
    else ofs << "-1\t";
    ofs << bb.upper(0) <<"\t";
    if (bb.rank>1) ofs << bb.upper(1) <<"\t";
    else ofs << "-1\t";
    if (bb.rank>2) ofs << bb.upper(2) <<"\t";
    else ofs << "-1\t";
    ofs << std::endl;
  }

protected:
  double last[3];
  int TraceDumpEvery;