00001
00002
00003
00004
00005
00006 #ifndef AMROC_CLP_FIXUP_H
00007 #define AMROC_CLP_FIXUP_H
00008
00016 #include "AMRFixup.h"
00017
00018
00019
00020 #include <iostream>
00021 #include <cstdlib>
00022
00023 template <class VectorType, class AuxVectorType, int dim> class ClpIntegrator;
00024
00031 template <class VectorType, int dim>
00032 class ClpFixupVecOps : private FixupOps<VectorType,VectorType,dim> {
00033 public:
00034 ClpFixupVecOps() {}
00035 protected:
00036 inline void copyv_to(GridData<VectorType,minus_1<dim>::dim> &target,
00037 const GridData<VectorType,dim> &source,
00038 const BBox &where, const int s)
00039 { FixupOps<VectorType,VectorType,dim>::copy_to(target, source, where, s); }
00040 };
00041
00042 #define f_rcflx FORTRAN_NAME(rcflx, RCFLX)
00043 extern "C" {
00044 void f_rcflx(const INTEGER& mflx);
00045 }
00046
00053 template <class VectorType, class FixupType, class AuxVectorType, int dim>
00054 class ClpFixup :
00055 public AMRFixup<VectorType,FixupType,dim>,
00056 public ClpFixupVecOps<VectorType,dim> {
00057
00058 typedef typename VectorType::InternalDataType DataType;
00059 typedef AMRFixup<VectorType,FixupType,dim> base;
00060 typedef ClpFixupVecOps<VectorType,dim> ops;
00061
00062 public:
00063 typedef ClpIntegrator<VectorType,AuxVectorType,dim> integrator_type;
00064 typedef typename base::vec_grid_data_type vec_grid_data_type;
00065 typedef typename base::ld_vec_grid_data_type ld_vec_grid_data_type;
00066
00067 ClpFixup(integrator_type& integ) : base(), _Integrator(integ) {}
00068
00069 protected:
00070 void CorrectFirstOrder(const int Time, const int Level, const int c,
00071 const double tc, const double tf,
00072 const double dtc, const double dtf,
00073 const int& mdim) {
00074
00075 int TimeC = CurrentTime(base::GH(),Level-1);
00076 DCoords dx = base::GH().worldStep(base::U().GF_StepSize(Level));
00077 const int refinefactor = RefineFactor(base::GH(),Level-1);
00078 #ifdef DEBUG_PRINT_WHOLE_FIXUP
00079 ( comm_service::log() << "\n*** Adding First Order Fluxes - "
00080 << base::U().interiorbbox(Time,Level,c) << " ***\n").flush();
00081 #endif
00082
00083 double Fac[2*dim];
00084 int d;
00085 for (d=0; d<dim; d++) {
00086 Fac[2*d] = dtf; Fac[2*d+1] = -dtf;
00087 for (int dd=0; dd<dim; dd++)
00088 if (d != dd) {
00089 Fac[2*d] *= dx(dd); Fac[2*d+1] *= dx(dd);
00090 }
00091 }
00092
00093 for (d=0; d<2*dim; d++)
00094 if ((mdim==0 || d/2==mdim-1) && base::ValidSide(Time,Level,c,d)) {
00095 #ifdef DEBUG_PRINT_WHOLE_FIXUP
00096 ( comm_service::log() << "Side: " << d << "\n").flush();
00097 #endif
00098 BBox bb(base::StateVecBBox(base::U().interiorbbox(Time,Level,c), d));
00099 BBox bbC(base::StateVecBBox(coarsen(base::U().interiorbbox(Time,Level,c),
00100 refinefactor),d));
00101 BBox workbb(bb); base::AlignBBox(workbb, d);
00102
00103 ld_vec_grid_data_type LeftState(workbb);
00104 ld_vec_grid_data_type RightState(workbb);
00105 ld_vec_grid_data_type NewState(workbb);
00106 BBox bbl(bb), bbr(bb);
00107 double dtl = dtf, dtr = dtf, tl, tr;
00108
00109
00110 if (d % 2 == 0) {
00111
00112 #ifdef DEBUG_AMR
00113 if (Integrator_().Abort())
00114 Integrator_().CheckGrid(base::U()(Time,Level,c), Level, bb, tf,
00115 "ClpFixup2::CorrectFirstOrder::LeftState,fine");
00116 #endif
00117 ops::copyv_to(LeftState, base::U()(Time,Level,c), bb, d);
00118 tl = tf;
00119
00120
00121 for (register int j=0; j<base::U().parents(Time,Level,c) ; j++) {
00122 BBox wbbC(coarsen(base::U().parentbox(Time,Level,c,j),refinefactor)*bbC);
00123 if (!wbbC.empty()) {
00124 #ifdef DEBUG_AMR
00125 if (Integrator_().Abort())
00126 Integrator_().CheckGrid(base::U()(TimeC,Level-1,base::U().parentidx(Time,Level,c,j)), Level-1,
00127 wbbC, tc, "ClpFixup2::CorrectFirstOrder::RightState,coarse");
00128 #endif
00129 ops::copyv_to(RightState, base::U()(TimeC,Level-1,base::U().parentidx(Time,Level,c,j)), wbbC, d);
00130 }
00131 }
00132 tr = tc;
00133 dtr = dtc;
00134 bbr.coarsen(d/2, refinefactor);
00135 }
00136 else {
00137
00138 for (register int j=0; j<base::U().parents(Time,Level,c) ; j++) {
00139 BBox wbbC(coarsen(base::U().parentbox(Time,Level,c,j),refinefactor)*bbC);
00140 if (!wbbC.empty()) {
00141 #ifdef DEBUG_AMR
00142 if (Integrator_().Abort())
00143 Integrator_().CheckGrid(base::U()(TimeC,Level-1,base::U().parentidx(Time,Level,c,j)), Level-1,
00144 wbbC, tc, "ClpFixup2::CorrectFirstOrder::LeftState,coarse");
00145 #endif
00146 ops::copyv_to(LeftState, base::U()(TimeC,Level-1,base::U().parentidx(Time,Level,c,j)), wbbC, d);
00147 }
00148 }
00149 tl = tc;
00150 dtl = dtc;
00151 bbl.coarsen(d/2, refinefactor);
00152
00153
00154 #ifdef DEBUG_AMR
00155 if (Integrator_().Abort())
00156 Integrator_().CheckGrid(base::U()(Time,Level,c), Level, bb, tf,
00157 "ClpFixup2::CorrectFirstOrder::RightState,fine");
00158 #endif
00159 ops::copyv_to(RightState, base::U()(Time,Level,c), bb, d);
00160 tr = tf;
00161 }
00162
00163 ((integrator_type&)Integrator_()).CalculateRP(LeftState, bbl, tl, dtl,
00164 RightState, bbr, tr, dtr, NewState, d/2+1);
00165
00166 #ifdef DEBUG_PRINT_WHOLE_FIXUP
00167 base::debug_print_ld(base::F(d)(Time,Level,c));
00168 ( comm_service::log() << "\n").flush();
00169 ( comm_service::log() << "LeftState: \n").flush();
00170 base::debug_print_ldv(LeftState); ( comm_service::log() << "\n").flush();
00171 ( comm_service::log() << "RightState: \n").flush();
00172 base::debug_print_ldv(RightState); ( comm_service::log() << "\n").flush();
00173 ( comm_service::log() << "NewState: \n").flush();
00174 base::debug_print_ldv(NewState); ( comm_service::log() << "\n").flush();
00175 #endif
00176 NewState.multiply(Fac[d],workbb);
00177 base::addf_to(base::F(d)(Time,Level,c),NewState,workbb);
00178 #ifdef DEBUG_PRINT_WHOLE_FIXUP
00179 base::debug_print_ld(base::F(d)(Time,Level,c));
00180 ( comm_service::log() << "\n").flush();
00181 #endif
00182 }
00183 }
00184
00185 virtual void AddFluxes(const int Time, const int Level, const int c,
00186 vec_grid_data_type* flux[],
00187 const double tc, const double tf,
00188 const double dtc, const double dtf,
00189 const int& mdim) {
00190 int mflx;
00191 f_rcflx(mflx);
00192 if (mflx==0) {
00193 if (dim>1 && Integrator_().DimSplit()) {
00194 std::cerr << "\n\nThe fixup is incorrect, if dimensional splitting is used with a\n";
00195 std::cerr << "scheme returning flux differences. Set Method(3)>=0 instead.\n";
00196 std::cerr << "Aborting programm...\n" << std::flush;
00197 std::exit(-1);
00198 }
00199 CorrectFirstOrder(Time, Level, c, tc, tf, dtc, dtf, mdim);
00200 }
00201
00202 base::AddIntercellFluxes(Time, Level, c, flux, dtf, mdim);
00203 }
00204
00205 inline integrator_type& Integrator_() { return _Integrator; }
00206 inline const integrator_type& Integrator_() const { return _Integrator; }
00207 inline const int& NGhosts() const { return _Integrator.NGhosts(); }
00208
00209 private:
00210 integrator_type& _Integrator;
00211 };
00212
00213 #endif