clawpack/applications/euler_znd/1d/StatDet/src/Problem.h

// -*- C++ -*-

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

#ifndef AMROC_PROBLEM_H
#define AMROC_PROBLEM_H

#include "eulerznd1.h"

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

class ExactSolutionSpecific : public F77ExactSolution<VectorType,DIM> {
  typedef VectorType::InternalDataType DataType;
  typedef F77ExactSolution<VectorType,DIM> base;
  typedef AMRSolver<VectorType,FixupType,FlagType,DIM> solver_type;
  typedef F77FileOutput<VectorType,DIM> output_type;
public:
  typedef base::vec_grid_fct_type vec_grid_fct_type;  
  typedef base::vec_grid_data_type vec_grid_data_type;
  typedef base::generic_func_type generic_func_type;

  ExactSolutionSpecific(solver_type& s, generic_func_type exact) : 
    base(exact), solver(s) {
    base::ENormAll = 1;
  }

  virtual void ErrorNorm(vec_grid_fct_type& u, grid_fct_type& work, 
                         const double& t) {
    if (ENorm<1 || ENorm>3) 
      return;
    int me = MY_PROC;
    for (int lev=0; lev<=FineLevel(base::GH()); lev++) {
      int Time = CurrentTime(base::GH(),lev); 
      int TStep = TimeStep(u,lev);
      forall (u,Time+TStep,lev,c) 
        SetGrid(u(Time+TStep,lev,c), lev, t);
      end_forall      
      int density_cnt = 1;
      ((output_type&) solver.FileOutput_()).Transform(u, work, Time, lev, density_cnt, t);
      ((output_type&) solver.FileOutput_()).Transform(u, work, Time+TStep, lev, density_cnt, t);
      forall (work,Time+TStep,lev,c) 
        work(Time+TStep,lev,c).minus(work(Time,lev,c));
      end_forall            
    }

    DataType *Error = (DataType *) 0;
    CalculateNorm(work, Error);

    if (me == VizServer) {
      std::ofstream outfile;
      std::ostream* out;
      std::string fname = "error.dat";
      outfile.open(fname.c_str(), std::ios::out | std::ios::app);
      out = new std::ostream(outfile.rdbuf());
      *out << t << " " << Error[0]<< std::endl;;
      outfile.close();
      delete out;
    }

    if (Error) delete [] Error;
  }  
protected:
  solver_type solver;
};
  

class SolverSpecific : 
  public AMRSolver<VectorType,FixupType,FlagType,DIM> {
  typedef VectorType::InternalDataType DataType;
  typedef AMRSolver<VectorType,FixupType,FlagType,DIM> base;
  typedef F77FileOutput<VectorType,DIM> output_type;
public:
  SolverSpecific(IntegratorSpecific& integ, 
                 base::initial_condition_type& init,
                 base::boundary_conditions_type& bc) :
  AMRSolver<VectorType,FixupType,FlagType,DIM>(integ, init, bc), OutputPMax(0) {
    SetLevelTransfer(new F77LevelTransfer<VectorType,DIM>(f_prolong, f_restrict));
    SetFileOutput(new F77FileOutput<VectorType,DIM>(f_flgout)); 
    SetFixup(new FixupSpecific(integ));
    SetFlagging(new FlaggingSpecific(*this)); 
    SetExactSolution(new ExactSolutionSpecific(*this,f_exact)); 
  }  

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

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

  virtual void Advance(double& t, double& dt) {
    base::Advance(t,dt);

    if (OutputPMax) {
      double maxp[2];
      for (register int l=0; l<=FineLevel(base::GH()); l++) {
        int Time = CurrentTime(base::GH(),l);
        int press_cnt = base::Dim()+4;
        ((output_type*) _FileOutput)->Transform(base::U(), base::Work(), Time, l, 
                                                press_cnt, base::t[l]);
        maxp[0] = 0.0; maxp[1] = 0.0; 
        forall (base::Work(),Time,l,c)       
          BBox bb = base::Work()(Time,l,c).bbox();
          BeginFastIndex1(work, bb, base::Work()(Time,l,c).data(), DataType);
          for (int n=bb.upper(0)-bb.stepsize(0); n>=bb.lower(0)+bb.stepsize(0); n-=bb.stepsize(0)) 
            if (FastIndex1(work,n)>FastIndex1(work,n+bb.stepsize(0))+1e-10 &&
                FastIndex1(work,n)>FastIndex1(work,n-bb.stepsize(0))+1e-10) {
              if (n>maxp[0]) { maxp[0] = double(n); maxp[1] = FastIndex1(work,n); }
              break;
            }
          EndFastIndex1(work);
        end_forall 
      }
#ifdef DAGH_NO_MPI
#else
      if (comm_service::dce() && comm_service::proc_num() > 1) {        
        int num = comm_service::proc_num();       
        int R;
        int size = 2*sizeof(double);
        void *values = (void *) new char[size*num];     
        R = MPI_Allgather(&maxp, size, MPI_BYTE, values, size, MPI_BYTE, 
                          comm_service::comm());
        if ( MPI_SUCCESS != R ) 
          comm_service::error_die( "SolverControlSpecific::maxp", "MPI_Allgather", R );
        for (int i=0;i<num;i++) {
          double tmaxn = *((double *)values + 2*i);
          double tmaxp = *((double *)values + 2*i+1);
          if (tmaxn>maxp[0]) maxp[1] = tmaxp;
        }
        if (values) delete [] (char *)values;
      }
#endif
      int me = MY_PROC; 
      if (me == VizServer) {
        std::ofstream outfile;
        std::ostream* out;
        std::string fname = "pmax.dat";
        outfile.open(fname.c_str(), std::ios::out | std::ios::app);
        out = new std::ostream(outfile.rdbuf());
        *out << t << " " << maxp[1] << std::endl;   
        outfile.close();
        delete out;
      }      
    }    
  }

protected:
  int OutputPMax;