00001
00002
00003
00004
00005
00006 #ifndef AMROC_AMRFIXUP_H
00007 #define AMROC_AMRFIXUP_H
00008
00016 #include "Fixup1.h"
00017 #include "Fixup2.h"
00018 #include "Fixup3.h"
00019
00020 #include "Timing.h"
00021
00022 #include <string>
00023
00044 template <class VectorType, class FixupType, int dim>
00045 class AMRFixup : public Fixup<VectorType,FixupType,dim>,
00046 public FixupOps<VectorType,FixupType,dim> {
00047
00048 typedef Fixup<VectorType,FixupType,dim> base;
00049 typedef FixupOps<VectorType,FixupType,dim> ops;
00050 typedef typename VectorType::InternalDataType DataType;
00051
00052 public:
00053 typedef GridData<VectorType,dim> vec_grid_data_type;
00054
00055 typedef GridData<VectorType,minus_1<dim>::dim> ld_vec_grid_data_type;
00056 typedef GridData<FixupType,minus_1<dim>::dim> ld_fixup_grid_data_type;
00057
00058 typedef GridFunction<VectorType,dim> vec_grid_fct_type;
00059 typedef GridFunction<FixupType,minus_1<dim>::dim> ld_fixup_grid_fct_type;
00060
00061 AMRFixup() : base() {
00062 _u = (vec_grid_fct_type *) 0;
00063 sprintf(GFName,"u-fixup");
00064 }
00065
00066 virtual ~AMRFixup() {}
00067
00068 virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00069 base::LocCtrl = Ctrl.getSubDevice(prefix+"AMRFixup");
00070 RegisterAt(base::LocCtrl,"Equations",base::_FixupEquations);
00071 }
00072 virtual void register_at(ControlDevice& Ctrl) { register_at(Ctrl, ""); }
00073
00074 virtual void SaveFluxes(const int Time, const int Level, const int c,
00075 vec_grid_data_type* flux[],
00076 const double dt, const int& mdim) {
00077
00078 if (!U().has_children(Time,Level,c)) return;
00079 DCoords dx = base::GH().worldStep(U().GF_StepSize(Level));
00080 const int refinefactor = RefineFactor(base::GH(),Level);
00081
00082 int d;
00083 DataType Fac = -dt;
00084 for (d=0; d<dim; d++)
00085 Fac *= dx(d);
00086
00087 #ifdef DEBUG_PRINT_WHOLE_FIXUP
00088 ( comm_service::log() << "\n*** Saving Coarse Fluxes - "
00089 << U().interiorbbox(Time,Level,c) << " ***\n").flush();
00090 #endif
00091 for (register int j=0; j<U().children(Time,Level,c); j++) {
00092 BBox work(U().childbox(Time,Level,c,j)*U().interiorbbox(Time,Level,c));
00093 if (!work.empty()) {
00094 const int& idx = U().childidx(Time,Level,c,j);
00095 for (d=0; d<2*dim; d++)
00096 if ((mdim==0 || d/2==mdim-1) && ValidSide(Time,Level+1,idx,d)) {
00097 BBox bb(FixupBBox(coarsen(U().interiorbbox(Time,Level+1,idx),refinefactor),d)*work);
00098 if (!bb.empty()) {
00099 bb = FluxBBox(bb,d);
00100 #ifdef DEBUG_PRINT_WHOLE_FIXUP
00101 ( comm_service::log() << "To " << U().interiorbbox(Time,Level+1,idx)
00102 << " Side: " << d << " using " << bb << " into "
00103 << base::F(d)(Time,Level+1,idx).bbox() << "\n").flush();
00104 #endif
00105 (*flux[d]).multiply(Fac/dx(d/2), bb);
00106 copyf_to(base::F(d)(Time,Level+1,idx),(*flux[d]),bb,d);
00107 (*flux[d]).divide(Fac/dx(d/2), bb);
00108 #ifdef DEBUG_PRINT_WHOLE_FIXUP
00109 ops::AlignBBox(bb,d);
00110 base::debug_print_ld(base::F(d)(Time,Level+1,idx),bb);
00111 ( comm_service::log() << "\n" ).flush();
00112 #endif
00113 }
00114 }
00115 }
00116 }
00117 }
00118
00119 virtual void AddFluxes(const int Time, const int Level, const int c,
00120 vec_grid_data_type* flux[],
00121 const double tc, const double tf,
00122 const double dt, const int& mdim) {
00123 AddIntercellFluxes(Time, Level, c, flux, dt, mdim);
00124 }
00125
00126 virtual void AddIntercellFluxes(const int Time, const int Level, const int c,
00127 vec_grid_data_type* flux[], const double dt,
00128 const int& mdim) {
00129
00130 DCoords dx = base::GH().worldStep(U().GF_StepSize(Level));
00131
00132 int d;
00133 DataType Fac = dt;
00134 for (d=0; d<dim; d++)
00135 Fac *= dx(d);
00136
00137 #ifdef DEBUG_PRINT_WHOLE_FIXUP
00138 ( comm_service::log() << "\n*** Adding Higher Order Fluxes - "
00139 << U().interiorbbox(Time,Level,c) << " ***\n").flush();
00140 #endif
00141 for (d=0; d<2*dim; d++)
00142 if ((mdim==0 || d/2==mdim-1) && ValidSide(Time,Level,c,d)) {
00143 BBox bb(FluxBBox(U().interiorbbox(Time,Level,c), d));
00144 (*flux[d]).multiply(Fac/dx(d/2), bb);
00145 addf_to(base::F(d)(Time,Level,c),(*flux[d]),bb,d);
00146 (*flux[d]).divide(Fac/dx(d/2), bb);
00147 #ifdef DEBUG_PRINT_WHOLE_FIXUP
00148 base::debug_print_ld(base::F(d)(Time,Level,c));
00149 ( comm_service::log() << "\n").flush();
00150 #endif
00151 }
00152 }
00153
00154 virtual void Correction(const int Time, const int WTime, const int Level) {
00155 int t_sten = 0;
00156 int s_sten = 1;
00157 vec_grid_fct_type uh(GFName, t_sten, s_sten, CurrentTime(base::GH(),Level), Level,
00158 base::GH(), DAGHCellCentered, 0, 1, DAGHCommSimple,
00159 DAGHNoBoundary, DAGHNoAdaptBoundary, DAGHNoExternalGhost);
00160
00161 #ifdef DEBUG_PRINT_WHOLE_FIXUP
00162 ( comm_service::log() << "\n*********** Conservative Fixup ***********\n").flush();
00163 #endif
00164
00165 for (int d=0; d<2*dim; d++)
00166 ParallelFixup(Time, Level, d, uh);
00167 }
00168
00169 inline void SetGridFunctions(vec_grid_fct_type* u) { _u = u; }
00170 inline vec_grid_fct_type& U() { return *_u; }
00171 inline const vec_grid_fct_type& U() const { return *_u; }
00172
00173 protected:
00174 int ValidSide(const int Time, const int Level, const int c, const int s) const {
00175 return (U().has_adaptiveboundary(Time,Level,c,s));
00176 }
00177
00178 void ParallelFixup(const int Time, const int Level,
00179 const int d, vec_grid_fct_type& uh) {
00180
00181 int TimeF = CurrentTime(base::GH(),Level+1);
00182 int TimeC = CurrentTime(base::GH(),Level);
00183 const int refinefactor = RefineFactor(base::GH(),Level);
00184 DCoords dx = base::GH().worldStep(U().GF_StepSize(Level));
00185
00186 DataType Fac = 1.0;
00187 for (int dd=0; dd<dim; dd++)
00188 Fac *= dx(dd);
00189
00190 forall (uh,TimeC,Level,c)
00191 uh(TimeC,Level,c) = 0.0;
00192 end_forall
00193
00194 forall (U(),TimeF,Level+1,c)
00195 if (ValidSide(TimeF,Level+1,c,d)) {
00196 BBox to(FixupBBox(coarsen(U().interiorbbox(TimeF,Level+1,c),refinefactor),d));
00197 for (register int j=0; j<U().parents(TimeF,Level+1,c) ; j++) {
00198 BBox wto(coarsen(U().parentbox(TimeF,Level+1,c,j),refinefactor)*to);
00199 const int& idx = U().parentidx(TimeF,Level+1,c,j);
00200 wto *= uh(TimeC,Level,idx).bbox();
00201 if (!wto.empty()) {
00202 copyf_from(uh(TimeC,Level,idx),wto,base::F(d)(TimeF,Level+1,c),d);
00203 uh(TimeC,Level,idx).multiply((d%2==0 ? -1.0 : 1.0) / Fac,wto);
00204 #ifdef DEBUG_PRINT_WHOLE_FIXUP
00205 ( comm_service::log() << "Writing at Side: " << d << " of "
00206 << U().interiorbbox(TimeF,Level+1,c) << " using " << wto << "\n").flush();
00207 base::debug_print(uh(TimeC,Level,idx), wto);
00208 ( comm_service::log() << "\n" ).flush();
00209 #endif
00210 }
00211 }
00212 }
00213 end_forall
00214
00215 END_WATCH_FIXUP_WHOLE
00216 START_WATCH
00217 Sync(uh,TimeC,Level);
00218 END_WATCH_FIXUP_SYNC
00219 START_WATCH
00220
00221 forall (U(),Time,Level,c)
00222 BBox to = U().interiorbbox(Time,Level,c);
00223 BBox from(to);
00224 from.shift(d/2, d%2 ? -1 : 1);
00225 U()(Time,Level,c).plus(uh(TimeC,Level,c), to, from);
00226 end_forall
00227 }
00228
00229 BBox StateVecBBox(const BBox& interior, const int s) {
00230 int d = s/2;
00231 #ifdef DEBUG_PRINT
00232 assert (d>=0 && d<dim);
00233 #endif
00234 BBox bb(interior);
00235 if (s%2 == 0) {
00236 bb.growlower(d,-1);
00237 bb.upper(d) = bb.lower(d);
00238 }
00239 else {
00240 bb.growupper(d,1);
00241 bb.lower(d) = bb.upper(d);
00242 }
00243 return bb;
00244 }
00245
00246 BBox FixupBBox(const BBox& interior, const int s) {
00247 int d = s/2;
00248 #ifdef DEBUG_PRINT
00249 assert (d>=0 && d<dim);
00250 #endif
00251 BBox bb(interior);
00252 if (s%2 == 0)
00253 bb.upper(d) = bb.lower(d);
00254 else
00255 bb.lower(d) = bb.upper(d);
00256 return bb;
00257 }
00258
00259 virtual BBox FluxBBox(const BBox& interior, const int s) {
00260 int d = s/2;
00261 #ifdef DEBUG_PRINT
00262 assert (d>=0 && d<dim);
00263 #endif
00264 BBox bb(interior);
00265 if (s%2 == 0)
00266 bb.upper(d) = bb.lower(d);
00267 else {
00268 bb.growupper(d,1);
00269 bb.lower(d) = bb.upper(d);
00270 }
00271 return bb;
00272 }
00273
00274 protected:
00275 inline void copyf_to(ld_fixup_grid_data_type &target, const BBox &to,
00276 const vec_grid_data_type &source,
00277 const BBox &from, const int s)
00278 { ops::copy_to(target, to, source, from, s); }
00279 inline void copyf_to(ld_fixup_grid_data_type &target,
00280 const vec_grid_data_type &source,
00281 const BBox &where, const int s)
00282 { ops::copy_to(target, source, where, s); }
00283
00284 inline void copyf_from(vec_grid_data_type &target, const BBox &to,
00285 const ld_fixup_grid_data_type &source,
00286 const BBox &from, const int s)
00287 { ops::copy_from(target, to, source, from, s); }
00288 inline void copyf_from(vec_grid_data_type &target, const BBox &where,
00289 const ld_fixup_grid_data_type &source, const int s)
00290 { ops::copy_from(target, where, source, s); }
00291
00292 inline void addf_to(ld_fixup_grid_data_type &target, const BBox &to,
00293 const vec_grid_data_type &source,
00294 const BBox &from, const int s)
00295 { ops::add_to(target, to, source, from, s); }
00296 inline void addf_to(ld_fixup_grid_data_type &target,
00297 const vec_grid_data_type &source,
00298 const BBox &where, const int s)
00299 { ops::add_to(target, source, where, s); }
00300
00301 inline void addf_from(vec_grid_data_type &target, const BBox &to,
00302 const ld_fixup_grid_data_type &source,
00303 const BBox &from, const int s)
00304 { ops::add_from(target, to, source, from, s); }
00305 inline void addf_from(vec_grid_data_type &target, const BBox &where,
00306 const ld_fixup_grid_data_type &source, const int s)
00307 { ops::add_from(target, where, source, s); }
00308
00309 inline void addf_to(ld_fixup_grid_data_type &target, const BBox &to,
00310 const ld_vec_grid_data_type &source,
00311 const BBox &from)
00312 { ops::add_to(target, to, source, from); }
00313 inline void addf_to(ld_fixup_grid_data_type &target,
00314 const ld_vec_grid_data_type &source,
00315 const BBox &where)
00316 { ops::add_to(target, source, where); }
00317
00318 protected:
00319 vec_grid_fct_type* _u;
00320 char GFName[DAGHBktGFNameWidth];
00321 };
00322
00323 #endif
