vtf-logo

clawpack/applications/euler_chem/1d/RealDetonation/src/Problem.h

00001 // -*- C++ -*-
00002 
00003 // Copyright (C) 2003-2007 California Institute of Technology
00004 // Ralf Deiterding, ralf@amroc.net
00005 
00006 #ifndef AMROC_PROBLEM_H
00007 #define AMROC_PROBLEM_H
00008 
00009 #include "eulerrhok1.h"
00010 #include "components.h"
00011 
00012 #define NAUX        0
00013 
00014 #include "ClpProblem.h"
00015 
00016 #define OWN_FLAGGING
00017 #define OWN_AMRSOLVER
00018 #include "ClpStdProblem.h"
00019 #include "F77Interfaces/F77UpdateLevelTransfer.h"
00020 
00021 class FlaggingSpecific : 
00022   public AMRFlagging<VectorType,FixupType,FlagType,DIM> {
00023   typedef AMRFlagging<VectorType,FixupType,FlagType,DIM> base;
00024 public:
00025   FlaggingSpecific(base::solver_type& solver) : base(solver) {
00026       base::AddCriterion(new ScaledGradient<VectorType,FlagType,DIM>());
00027       base::AddCriterion(new LimiterType<VectorType,FlagType,DIM>());
00028       base::AddCriterion(new AbsoluteError<VectorType,FixupType,FlagType,DIM>(solver));
00029       base::AddCriterion(new RelativeError<VectorType,FixupType,FlagType,DIM>(solver));
00030       base::AddCriterion(new F77ScaledGradient<VectorType,FlagType,DIM>(f_flg));
00031       base::AddCriterion(new F77LimiterType<VectorType,FlagType,DIM>(f_flg));
00032       base::AddCriterion(new F77AbsoluteError<VectorType,FixupType,FlagType,DIM>(solver,f_flg));
00033       base::AddCriterion(new F77RelativeError<VectorType,FixupType,FlagType,DIM>(solver,f_flg));
00034     }
00035 
00036   ~FlaggingSpecific() { DeleteAllCriterions(); }
00037 };
00038 
00039 
00040 class SolverSpecific : 
00041   public AMRSolver<VectorType,FixupType,FlagType,DIM> {
00042   typedef VectorType::InternalDataType DataType;
00043   typedef AMRSolver<VectorType,FixupType,FlagType,DIM> base;
00044   typedef F77FileOutput<VectorType,DIM> output_type;
00045 public:
00046   SolverSpecific(IntegratorSpecific& integ, 
00047                  base::initial_condition_type& init,
00048                  base::boundary_conditions_type& bc) :
00049     AMRSolver<VectorType,FixupType,FlagType,DIM>(integ, init, bc), OutputPMax(0), PThreshold(0.) {
00050       SetLevelTransfer(new F77UpdateLevelTransfer<VectorType,DIM>(f_prolong, f_restrict, f_tupdate));
00051     SetFileOutput(new output_type(f_out)); 
00052     SetFixup(new FixupSpecific(integ));
00053     SetFlagging(new FlaggingSpecific(*this)); 
00054   }  
00055 
00056   ~SolverSpecific() {
00057     delete _LevelTransfer;
00058     delete _Flagging;
00059     delete _Fixup;
00060     delete _FileOutput;
00061   }
00062 
00063   virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00064     base::register_at(Ctrl, prefix);
00065     RegisterAt(LocCtrl,"OutputPMax",OutputPMax); 
00066     RegisterAt(LocCtrl,"PThreshold",PThreshold); 
00067   }
00068   virtual void register_at(ControlDevice& Ctrl) {
00069     register_at(Ctrl, "");
00070   }
00071 
00072   virtual void AfterLevelStep(const int Level) {
00073     if (Level<FineLevel(base::GH())) return;
00074     if (OutputPMax) {
00075       double maxp[2];
00076       for (register int l=0; l<=FineLevel(base::GH()); l++) {
00077         int Time = CurrentTime(base::GH(),l);
00078         int press_cnt = base::Dim()+4;
00079         ((output_type*) _FileOutput)->Transform(base::U(), base::Work(), Time, l, 
00080                                                 press_cnt, base::t[l]);
00081         maxp[0] = 0.0; maxp[1] = 0.0; 
00082         forall (base::Work(),Time,l,c)       
00083           BBox bb = base::Work()(Time,l,c).bbox();
00084           BeginFastIndex1(work, bb, base::Work()(Time,l,c).data(), DataType);
00085           for (int n=bb.upper(0)-bb.stepsize(0); n>=bb.lower(0)+bb.stepsize(0); n-=bb.stepsize(0)) 
00086             if (FastIndex1(work,n)>FastIndex1(work,n+bb.stepsize(0)) &&
00087                 FastIndex1(work,n)>FastIndex1(work,n-bb.stepsize(0)) &&
00088                 FastIndex1(work,n)>PThreshold) {
00089               if (n>maxp[0]) { maxp[0] = double(n); maxp[1] = FastIndex1(work,n); }
00090               break;
00091             }
00092           EndFastIndex1(work);
00093         end_forall 
00094       }
00095 #ifdef DAGH_NO_MPI
00096 #else
00097       if (comm_service::dce() && comm_service::proc_num() > 1) {        
00098         int num = comm_service::proc_num();       
00099         int R;
00100         int size = 2*sizeof(double);
00101         void *values = (void *) new char[size*num];     
00102         R = MPI_Allgather(&maxp, size, MPI_BYTE, values, size, MPI_BYTE, 
00103                           comm_service::comm());
00104         if ( MPI_SUCCESS != R ) 
00105           comm_service::error_die( "SolverControlSpecific::maxp", "MPI_Allgather", R );
00106         for (int i=0;i<num;i++) {
00107           double tmaxn = *((double *)values + 2*i);
00108           double tmaxp = *((double *)values + 2*i+1);
00109           if (tmaxn>maxp[0]) maxp[1] = tmaxp;
00110         }
00111         if (values) delete [] (char *)values;
00112       }
00113 #endif
00114       int me = MY_PROC; 
00115       if (me == VizServer) {
00116         std::ofstream outfile;
00117         std::ostream* out;
00118         std::string fname = "pmax.dat";
00119         outfile.open(fname.c_str(), std::ios::out | std::ios::app);
00120         out = new std::ostream(outfile.rdbuf());
00121         *out << base::t[Level] << " " << maxp[1] << std::endl;   
00122         outfile.close();
00123         delete out;
00124       }      
00125     }    
00126   }
00127 
00128 protected:
00129   int OutputPMax;
00130   double PThreshold;
00131 };