00001
00002
00003
00004
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 };