00001
00002
00003
00004
00005
00006 #ifndef AMROC_PROBLEM_H
00007 #define AMROC_PROBLEM_H
00008
00009 #include "eulerznd1.h"
00010
00011 #include "ClpProblem.h"
00012 #define OWN_AMRSOLVER
00013 #include "ClpStdProblem.h"
00014
00015 class ExactSolutionSpecific : public F77ExactSolution<VectorType,DIM> {
00016 typedef VectorType::InternalDataType DataType;
00017 typedef F77ExactSolution<VectorType,DIM> base;
00018 typedef AMRSolver<VectorType,FixupType,FlagType,DIM> solver_type;
00019 typedef F77FileOutput<VectorType,DIM> output_type;
00020 public:
00021 typedef base::vec_grid_fct_type vec_grid_fct_type;
00022 typedef base::vec_grid_data_type vec_grid_data_type;
00023 typedef base::generic_func_type generic_func_type;
00024
00025 ExactSolutionSpecific(solver_type& s, generic_func_type exact) :
00026 base(exact), solver(s) {
00027 base::ENormAll = 1;
00028 }
00029
00030 virtual void ErrorNorm(vec_grid_fct_type& u, grid_fct_type& work,
00031 const double& t) {
00032 if (ENorm<1 || ENorm>3)
00033 return;
00034 int me = MY_PROC;
00035 for (int lev=0; lev<=FineLevel(base::GH()); lev++) {
00036 int Time = CurrentTime(base::GH(),lev);
00037 int TStep = TimeStep(u,lev);
00038 forall (u,Time+TStep,lev,c)
00039 SetGrid(u(Time+TStep,lev,c), lev, t);
00040 end_forall
00041 int density_cnt = 1;
00042 ((output_type&) solver.FileOutput_()).Transform(u, work, Time, lev, density_cnt, t);
00043 ((output_type&) solver.FileOutput_()).Transform(u, work, Time+TStep, lev, density_cnt, t);
00044 forall (work,Time+TStep,lev,c)
00045 work(Time+TStep,lev,c).minus(work(Time,lev,c));
00046 end_forall
00047 }
00048
00049 DataType *Error = (DataType *) 0;
00050 CalculateNorm(work, Error);
00051
00052 if (me == VizServer) {
00053 std::ofstream outfile;
00054 std::ostream* out;
00055 std::string fname = "error.dat";
00056 outfile.open(fname.c_str(), std::ios::out | std::ios::app);
00057 out = new std::ostream(outfile.rdbuf());
00058 *out << t << " " << Error[0]<< std::endl;;
00059 outfile.close();
00060 delete out;
00061 }
00062
00063 if (Error) delete [] Error;
00064 }
00065 protected:
00066 solver_type solver;
00067 };
00068
00069
00070 class SolverSpecific :
00071 public AMRSolver<VectorType,FixupType,FlagType,DIM> {
00072 typedef VectorType::InternalDataType DataType;
00073 typedef AMRSolver<VectorType,FixupType,FlagType,DIM> base;
00074 typedef F77FileOutput<VectorType,DIM> output_type;
00075 public:
00076 SolverSpecific(IntegratorSpecific& integ,
00077 base::initial_condition_type& init,
00078 base::boundary_conditions_type& bc) :
00079 AMRSolver<VectorType,FixupType,FlagType,DIM>(integ, init, bc), OutputPMax(0) {
00080 SetLevelTransfer(new F77LevelTransfer<VectorType,DIM>(f_prolong, f_restrict));
00081 SetFileOutput(new F77FileOutput<VectorType,DIM>(f_flgout));
00082 SetFixup(new FixupSpecific(integ));
00083 SetFlagging(new FlaggingSpecific(*this));
00084 SetExactSolution(new ExactSolutionSpecific(*this,f_exact));
00085 }
00086
00087 ~SolverSpecific() {
00088 delete _LevelTransfer;
00089 delete _Flagging;
00090 delete _Fixup;
00091 delete _FileOutput;
00092 delete _ExactSolution;
00093 }
00094
00095 virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00096 base::register_at(Ctrl, prefix);
00097 RegisterAt(LocCtrl,"OutputPMax",OutputPMax);
00098 }
00099 virtual void register_at(ControlDevice& Ctrl) {
00100 register_at(Ctrl, "");
00101 }
00102
00103 virtual void Advance(double& t, double& dt) {
00104 base::Advance(t,dt);
00105
00106 if (OutputPMax) {
00107 double maxp[2];
00108 for (register int l=0; l<=FineLevel(base::GH()); l++) {
00109 int Time = CurrentTime(base::GH(),l);
00110 int press_cnt = base::Dim()+4;
00111 ((output_type*) _FileOutput)->Transform(base::U(), base::Work(), Time, l,
00112 press_cnt, base::t[l]);
00113 maxp[0] = 0.0; maxp[1] = 0.0;
00114 forall (base::Work(),Time,l,c)
00115 BBox bb = base::Work()(Time,l,c).bbox();
00116 BeginFastIndex1(work, bb, base::Work()(Time,l,c).data(), DataType);
00117 for (int n=bb.upper(0)-bb.stepsize(0); n>=bb.lower(0)+bb.stepsize(0); n-=bb.stepsize(0))
00118 if (FastIndex1(work,n)>FastIndex1(work,n+bb.stepsize(0))+1e-10 &&
00119 FastIndex1(work,n)>FastIndex1(work,n-bb.stepsize(0))+1e-10) {
00120 if (n>maxp[0]) { maxp[0] = double(n); maxp[1] = FastIndex1(work,n); }
00121 break;
00122 }
00123 EndFastIndex1(work);
00124 end_forall
00125 }
00126 #ifdef DAGH_NO_MPI
00127 #else
00128 if (comm_service::dce() && comm_service::proc_num() > 1) {
00129 int num = comm_service::proc_num();
00130 int R;
00131 int size = 2*sizeof(double);
00132 void *values = (void *) new char[size*num];
00133 R = MPI_Allgather(&maxp, size, MPI_BYTE, values, size, MPI_BYTE,
00134 comm_service::comm());
00135 if ( MPI_SUCCESS != R )
00136 comm_service::error_die( "SolverControlSpecific::maxp", "MPI_Allgather", R );
00137 for (int i=0;i<num;i++) {
00138 double tmaxn = *((double *)values + 2*i);
00139 double tmaxp = *((double *)values + 2*i+1);
00140 if (tmaxn>maxp[0]) maxp[1] = tmaxp;
00141 }
00142 if (values) delete [] (char *)values;
00143 }
00144 #endif
00145 int me = MY_PROC;
00146 if (me == VizServer) {
00147 std::ofstream outfile;
00148 std::ostream* out;
00149 std::string fname = "pmax.dat";
00150 outfile.open(fname.c_str(), std::ios::out | std::ios::app);
00151 out = new std::ostream(outfile.rdbuf());
00152 *out << t << " " << maxp[1] << std::endl;
00153 outfile.close();
00154 delete out;
00155 }
00156 }
00157 }
00158
00159 protected:
00160 int OutputPMax;