00001
00002
00003
00004
00005
00006
00007
00008
00009 #ifndef AMROC_AMRSOLVER_H
00010 #define AMROC_AMRSOLVER_H
00011
00019 #include "AMRSolverBase.h"
00020 #include "DAGHCluster.h"
00021
00022 #include <iostream>
00023 #include <sstream>
00024 #include <string>
00025 #include <cstdio>
00026 #include <cstdlib>
00027
00028 #define AMRSolverRegridFalse (0)
00029 #define VERY_LARGE_CFL (1.5)
00030 #define MAXOVERLAPWIDTH (100)
00031
00038
00039
00040
00041
00042
00043
00044
00045
00046 template <class VectorType, class FixupType, class FlagType, int dim>
00047 class AMRSolver : public AMRSolverBase<VectorType,FixupType,FlagType,dim> {
00048 typedef AMRSolverBase<VectorType,FixupType,FlagType,dim> base;
00049
00050 typedef Vector<FlagType,1> VFlagType;
00051 typedef GridData<VFlagType,dim> vflag_grid_data_type;
00052
00053 public:
00054 typedef GridData<FlagType,dim> flag_data_type;
00055 typedef typename base::vec_grid_fct_type vec_grid_fct_type;
00056 typedef typename base::vec_grid_data_type vec_grid_data_type;
00057 typedef typename base::grid_fct_type grid_fct_type;
00058 typedef typename base::grid_data_type grid_data_type;
00059
00060 typedef typename base::integrator_type integrator_type;
00061 typedef typename base::initial_condition_type initial_condition_type;
00062 typedef typename base::boundary_conditions_type boundary_conditions_type;
00063 typedef typename base::leveltransfer_type leveltransfer_type;
00064 typedef typename base::flagging_type flagging_type;
00065 typedef typename base::fixup_type fixup_type;
00066
00067 AMRSolver(integrator_type& integ, initial_condition_type& init,
00068 boundary_conditions_type& bc) : base(integ, init, bc),
00069 FixupPar(1), RegridEvery(1), MinEfficiency(0.7),
00070 BlockWidth(4), OverlapWidth(0), BufferWidth(2), NestingBuffer(1), NoTimeRefine(0),
00071 AdaptBndTimeInterpolate(1), PlotFlags(0) {
00072
00073 t = (double *)NULL; dt = (double *)NULL; cfl_new = (double *)NULL;
00074 AllocError::SetTexts("AMRSolver","Constructor");
00075 }
00076
00077 virtual ~AMRSolver() {
00078 if (t) delete [] t;
00079 if (dt) delete [] dt;
00080 if (cfl_new) delete [] cfl_new;
00081 }
00082
00083 virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00084 base::LocCtrl = Ctrl.getSubDevice(prefix+"AMRSolver");
00085 RegisterAt(base::LocCtrl,"DoFixup",FixupPar);
00086 RegisterAt(base::LocCtrl,"NoTimeRefine",NoTimeRefine);
00087 RegisterAt(base::LocCtrl,"AdaptBndTimeInterpolate",AdaptBndTimeInterpolate);
00088 RegisterAt(base::LocCtrl,"RegridEvery",RegridEvery);
00089 RegisterAt(base::LocCtrl,"MinEfficiency",MinEfficiency);
00090 RegisterAt(base::LocCtrl,"BufferWidth",BufferWidth);
00091 RegisterAt(base::LocCtrl,"BlockWidth",BlockWidth);
00092 RegisterAt(base::LocCtrl,"OverlapWidth",OverlapWidth);
00093 RegisterAt(base::LocCtrl,"NestingBuffer",NestingBuffer);
00094 RegisterAt(base::LocCtrl,"OutputFlags",PlotFlags);
00095 base::register_at(base::LocCtrl,"");
00096 }
00097 virtual void register_at(ControlDevice& Ctrl) {
00098 register_at(Ctrl, "");
00099 }
00100
00101 virtual void SetupData() {
00102 base::SetupData();
00103
00104 if (base::_Flagging && MaxLevel(base::GH())>=1)
00105 base::Flagging_().SetBufferwidth(std::min(std::max(BufferWidth,OverlapWidth-1),
00106 MAXOVERLAPWIDTH));
00107 else
00108 RegridEvery = 0;
00109 if (FixupPar == 0)
00110 base::SetFixup((fixup_type*) 0);
00111
00112 t = new double[MaxLevel(base::GH())+1];
00113 dt = new double[MaxLevel(base::GH())+1];
00114 cfl_new = new double[MaxLevel(base::GH())+1];
00115
00116 if (FixupPar && NestingBuffer<1 && RegridEvery!=AMRSolverRegridFalse) {
00117 if (MY_PROC == VizServer)
00118 std::cout << " Fixup requires a NestingBuffer>0 ! Using Default." << std::endl;
00119 NestingBuffer = 1;
00120 }
00121 base::GH().DAGH_SetNestingBuffer(NestingBuffer);
00122
00123 if (base::_Flagging)
00124 base::_Flagging->SetGridFunctions(base::_u, base::_u_sh, base::_work, base::_work_sh);
00125 if (base::_Fixup)
00126 base::_Fixup->SetGridFunctions(base::_u);
00127
00128 if (NoTimeRefine==1) SetTimeRefineWL(base::GH(),DAGHTrue);
00129 else SetTimeRefineWL(base::GH(),DAGHFalse);
00130 }
00131
00132 virtual void Initialize_(const double& dt_start) {
00133 START_WATCH
00134
00135 t[0] = 0.0; dt[0] = dt_start; cfl_new[0] = 0.0;
00136 for (register int l=1; l<=MaxLevel(base::GH()); l++) {
00137 t[l] = 0.0;
00138 cfl_new[l] = 0.0;
00139 }
00140
00141 base::U().GF_DeleteStorage(1);
00142 if (base::ErrorEstimation())
00143 base::Ush().GF_DeleteStorage(1);
00144
00145
00146
00147
00148 int lev = 0;
00149 if (RegridEvery != AMRSolverRegridFalse) {
00150 while (lev<=FineLevel(base::GH()) && lev<MaxLevel(base::GH())) {
00151
00152
00153
00154 #ifdef DEBUG_PRINT_AMRSOLVER
00155 ( comm_service::log() << " *** Setting flags at Level: " << lev << "\n" ).flush();
00156 #endif
00157
00158
00159
00160 for (register int l=0; l<=FineLevel(base::GH()); l++) {
00161 base::InitialCondition_().Set(base::U(), base::Work(), l);
00162 SetBndry(base::U(),0,lev,0.0);
00163 if (base::ErrorEstimation()) {
00164 base::InitialCondition_().Set(base::Ush(), base::Worksh(), l);
00165 SetBndry(base::Ush(),0,lev,0.0);
00166 }
00167 }
00168 BBoxList bblist;
00169 RegridLevel(0, 0, bblist, false);
00170
00171 END_WATCH_INITIALIZATION
00172 START_WATCH
00173 RecomposeHierarchy(base::GH());
00174 END_WATCH_RECOMPOSING_WHOLE
00175 START_WATCH
00176 lev++;
00177 }
00178
00179 if (FineLevel(base::GH())>0) {
00180 for (register int l=0; l<=FineLevel(base::GH()); l++) {
00181 base::InitialCondition_().Set(base::U(), base::Work(), l);
00182 SetBndry(base::U(),0,l,0.0);
00183 if (base::ErrorEstimation()) {
00184 base::InitialCondition_().Set(base::Ush(), base::Worksh(), l);
00185 SetBndry(base::Ush(),0,l,0.0);
00186 }
00187 }
00188 BBoxList bblist;
00189 forall (base::U(),0,FineLevel(base::GH()),c)
00190 bblist.add(coarsen(base::U().interiorbbox(0,FineLevel(base::GH()),c),
00191 RefineFactor(base::GH(),FineLevel(base::GH())-1)));
00192 end_forall
00193 RegridLevel(0, 0, bblist, true);
00194 END_WATCH_INITIALIZATION
00195 START_WATCH
00196 RecomposeHierarchy(base::GH());
00197 END_WATCH_RECOMPOSING_WHOLE
00198 START_WATCH
00199 }
00200 }
00201
00202 base::U().GF_CreateStorage(1);
00203 if (base::ErrorEstimation())
00204 base::Ush().GF_CreateStorage(1);
00205
00206 for (lev=0; lev<=FineLevel(base::GH()); lev++) {
00207 #ifdef DEBUG_PRINT_AMRSOLVER
00208 ( comm_service::log() << " *** Initializing Level: " << lev << " *** \n" ).flush();
00209 #endif
00210 t[lev] = 0.0;
00211 base::InitialCondition_().Set(base::U(), base::Work(), lev);
00212 SetBndry(base::U(),0,lev,t[lev]);
00213 if (base::ErrorEstimation()) {
00214 base::InitialCondition_().Set(base::Ush(), base::Worksh(), lev);
00215 SetBndry(base::Ush(),0,lev,t[lev]);
00216 }
00217 }
00218
00219 SetTimeSpecs(base::GH(), 0.0, dt_start, 2);
00220
00221 for (lev=FineLevel(base::GH()); lev>=1; lev--) {
00222 Restrict(base::U(), 0, lev, 0, lev-1);
00223 if (base::ErrorEstimation())
00224 Restrict(base::Ush(), 0, lev, 0, lev-1);
00225 }
00226
00227 END_WATCH_INITIALIZATION
00228 }
00229
00230 virtual bool Restart_(const char* CheckpointFile) {
00231 START_WATCH
00232 bool ret = RecomposeHierarchy(base::GH(), CheckpointFile);
00233 if (ret) {
00234 double NextTime;
00235 for (register int l=0; l<=FineLevel(base::GH()); l++) {
00236 base::GH().DAGH_GetTimeSpecs(t[l], NextTime);
00237 int Time = CurrentTime(base::GH(),l);
00238 SetPhysicalTime(base::U(),Time,l,t[l]);
00239 if (base::ErrorEstimation()) SetPhysicalTime(base::Ush(),Time,l,t[l]);
00240 }
00241 }
00242 END_WATCH_INITIALIZATION
00243 return ret;
00244 }
00245
00246 virtual void Checkpointing_(const char* CheckpointFile) {
00247 Checkpoint(base::GH(), CheckpointFile);
00248 }
00249
00250 virtual void Restart_(std::stringstream& CheckpointStr) {
00251 base::GH().DAGH_RecomposeHierarchy(CheckpointStr);
00252 }
00253
00254 virtual void Checkpointing_(std::stringstream& CheckpointStr) {
00255 base::GH().DAGH_Checkpoint(CheckpointStr);
00256 }
00257
00258 virtual double Tick(int VariableTimeStepping, const double dtv[],
00259 const double cflv[], int& Rejections) {
00260
00261 double told = t[0];
00262 double cfl_max;
00263 Rejections = 0;
00264 do {
00265 std::stringstream* RestartStorage = (std::stringstream *)NULL;
00266 if (VariableTimeStepping==AutomaticTimeStepsRestart) {
00267 AllocError::SetTexts("AMRSolver","Memory for restart");
00268 RestartStorage = new std::stringstream;
00269 Checkpointing_(*RestartStorage);
00270 }
00271
00272 #ifdef DEBUG_PRINT_AMRSOLVER
00273 double Current_Time, Next_Time;
00274 base::GH().DAGH_GetTimeSpecs(Current_Time, Next_Time);
00275 ( comm_service::log() << " --- Iteration: " << StepsTaken(base::GH(),0)+1 << " ---"
00276 << " t = " << Current_Time << " dt = "
00277 << DeltaT(base::GH(), 0) << "\n" ).flush();
00278 #endif
00279 #ifdef DEBUG_PROFILE
00280 ( std::cout << "\n" ).flush();
00281 #endif
00282
00283 for (register int l=1; l<=MaxLevel(base::GH()); l++)
00284 t[l] = t[0];
00285
00286 AdvanceLevel(0, RegridEvery,
00287 !(RegridEvery != AMRSolverRegridFalse ?
00288 (StepsTaken(base::GH(),0)%RegridEvery==0 && StepsTaken(base::GH(),0)>0) : true),
00289 base::ErrorEstimation() && (RegridEvery != AMRSolverRegridFalse ?
00290 (StepsTaken(base::GH(),0)+1)%RegridEvery==0 : false), FixupPar > 0,
00291 !(base::RedistributeEvery>0 && StepsTaken(base::GH(),0)%base::RedistributeEvery!=0),
00292 base::RedistributeEvery<=0);
00293
00294 cfl_max = 0.0;
00295 for (register int l=0; l<=FineLevel(base::GH()); l++)
00296 cfl_max = cfl_max > cfl_new[l] ? cfl_max : cfl_new[l];
00297
00298 #ifdef DAGH_NO_MPI
00299 #else
00300 if (comm_service::dce() && comm_service::proc_num() > 1 &&
00301 VariableTimeStepping >= 0) {
00302
00303 int num = comm_service::proc_num();
00304
00305 int R;
00306 int size = sizeof(double);
00307 void *values = (void *) new char[size*num];
00308
00309 START_WATCH
00310 R = MPI_Allgather(&cfl_max, size, MPI_BYTE, values, size, MPI_BYTE,
00311 comm_service::comm(base::U().GF_Id()));
00312 END_WATCH(GATHERCFL);
00313 if ( MPI_SUCCESS != R )
00314 comm_service::error_die( "AMRSolver::cfl_max", "MPI_Allgather", R );
00315
00316 for (int i=0;i<num;i++) {
00317 double tmax = *((double *)values + i);
00318 cfl_max = cfl_max > tmax ? cfl_max : tmax;
00319 }
00320 if (values) delete [] (char *)values;
00321 }
00322 #endif
00323
00324 #ifdef DEBUG_PRINT_AMRSOLVER
00325 ( comm_service::log() << " Levels = " << FineLevel(base::GH())+1
00326 << " max. CFL = " << cfl_max << "\n" ).flush();
00327 #endif
00328 #ifdef DEBUG_PROFILE
00329 ( std::cout << "\n" ).flush();
00330 #endif
00331
00332 if (cfl_max > cflv[0]) {
00333 if (VariableTimeStepping==AutomaticTimeStepsRestart) {
00334 if (RestartStorage!=(std::stringstream *)NULL) {
00335 Restart_(*RestartStorage);
00336 for (int lev=0; lev<=FineLevel(base::GH()); lev++)
00337 t[lev] = told;
00338 Rejections++;
00339 }
00340 else
00341 std::cerr << " Restart failure!" << std::endl << std::flush ;
00342 if (cfl_max > VERY_LARGE_CFL)
00343 std::cerr << " CFL>" << (double)VERY_LARGE_CFL << " " << std::flush;
00344 }
00345 else if (cfl_max > 1.0)
00346 std::cerr << " CFL>1.0 " << std::flush;
00347 }
00348
00349 if (VariableTimeStepping!=FixedTimeSteps && VariableTimeStepping>=0)
00350 if (cfl_max == 0.0)
00351 SetTimeSpecs(base::GH(), t[0], t[0]+dtv[1], 2);
00352 else {
00353 double dt_predict = (DeltaT(base::GH(), 0) / cfl_max) * cflv[1];
00354 SetTimeSpecs(base::GH(), t[0], t[0] +
00355 (dtv[1] < dt_predict ? dtv[1] : dt_predict), 2);
00356 }
00357 else
00358 SetTimeSpecs(base::GH(), t[0], t[0]+dtv[0], 2);
00359
00360 if (RestartStorage != (std::stringstream *)NULL)
00361 delete RestartStorage;
00362
00363 } while (t[0] <= told);
00364
00365 return cfl_max;
00366 }
00367
00368 virtual double IntegrateLevel(vec_grid_fct_type& u, const int Time, const int Level,
00369 double t, double dt, bool DoFixup, double tc,
00370 const int which) {
00371
00372 double cfl = 0, cfl_step;
00373 int mpass = 0;
00374 int TStep = TimeStep(u,Level);
00375 if (base::Integrator_().MaxIntegratorPasses() > 0) mpass = 1;
00376 do {
00377
00378 cfl_step = UpdateLevel(u,Time,Level,t,dt,DoFixup,tc,which,mpass);
00379 cfl = (cfl > cfl_step ? cfl : cfl_step);
00380
00381
00382
00383 if (mpass>0 && mpass<base::Integrator_().MaxIntegratorPasses()) {
00384 SwapTimeLevels(u,Level,Time,Time+TStep);
00385 SetBndry(u,Time,Level,t);
00386 SwapTimeLevels(u,Level,Time,Time+TStep);
00387 }
00388
00389 mpass++;
00390 } while (mpass <= base::Integrator_().MaxIntegratorPasses());
00391
00392 #ifdef DEBUG_PRINT_AMRSOLVER
00393 ( comm_service::log() << "IntegrateLevel(): Level" << Level << " local CFL = " << cfl << "\n" ).flush();
00394 #endif
00395 return cfl;
00396 }
00397
00398
00399
00400
00401
00402 virtual double UpdateLevel(vec_grid_fct_type& u, const int& Time, const int& Level,
00403 double t, double dtl, bool DoFixup, double tc,
00404 const int& which, const int& mpass) {
00405
00406 char errtext[256], whichtext[16];
00407 if (which<=0 || which>2) std::sprintf(whichtext,"Main");
00408 else if (which==1) std::sprintf(whichtext,"Estimate");
00409 else std::sprintf(whichtext,"Shadow");
00410 std::sprintf(errtext,"Flux allocation for %s",whichtext);
00411 AllocError::SetTexts("AMRSolver<dim>",errtext);
00412 base::Integrator_().SetGridFunction(&u);
00413
00414 int TStep = TimeStep(u,Level);
00415 double cfl = 0, cfl_grid;
00416
00417 forall (u,Time,Level,c)
00418 vec_grid_data_type& u_old = u(Time, Level, c);
00419 vec_grid_data_type& u_new = u(Time+TStep, Level, c);
00420
00421 vec_grid_data_type** f;
00422 base::Integrator_().AllocGridFluxes(u_old.bbox(), f);
00423
00424 #ifdef DEBUG_PRINT_INTEGRATOR
00425 if (mpass<=1) ( comm_service::log() << " " << which << ": " ).flush();
00426 #endif
00427 START_WATCH
00428 cfl_grid = base::Integrator_().CalculateGrid(u_new, u_old, f, Level, t, dtl, mpass);
00429 END_WATCH_INTEGRATION(which)
00430 #ifdef DEBUG_PROFILE
00431 if (mpass<=1) ( std::cout << "." ).flush();
00432 #endif
00433
00434 cfl = (cfl > cfl_grid ? cfl : cfl_grid);
00435
00436 if (Level>0 && DoFixup && cfl_grid>0 && cfl_grid<=VERY_LARGE_CFL) {
00437
00438
00439
00440 START_WATCH
00441 base::Fixup_().AddFluxes(Time, Level, c, f, tc, t, dt[Level-1], dtl, mpass);
00442 END_WATCH(FIXUP_WHOLE)
00443 }
00444
00445 if (Level<FineLevel(base::GH()) && DoFixup) {
00446 START_WATCH
00447
00448 if (!(cfl_grid>0 && cfl_grid<=VERY_LARGE_CFL))
00449 base::Integrator_().ResetGridFluxes(f);
00450 base::Fixup_().SaveFluxes(Time, Level, c, f , dtl, mpass);
00451 END_WATCH(FIXUP_WHOLE)
00452 }
00453
00454 base::Integrator_().DeAllocGridFluxes(f);
00455 end_forall
00456
00457 #ifdef DEBUG_PRINT_AMRSOLVER
00458 ( comm_service::log() << "UpdateLevel(): Level" << Level
00459 << " mpass = " << mpass << " local CFL = " << cfl << "\n" ).flush();
00460 #endif
00461 return cfl;
00462 }
00463
00464
00465
00466
00467
00468
00469
00470
00471
00472
00473
00474
00475
00476
00477
00478
00479
00480
00481
00482
00483
00484
00485
00486
00487 virtual void SetBndry(vec_grid_fct_type& u, const int Time,
00488 const int Level, double t) {
00489 START_WATCH
00490 AdaptiveBoundaryUpdate(u,Time,Level);
00491 END_WATCH(BOUNDARIES_INTERPOLATION)
00492 START_WATCH
00493 Sync(u,Time,Level);
00494 END_WATCH(BOUNDARIES_SYNC)
00495 START_WATCH
00496 ExternalBoundaryUpdate(u,Time,Level);
00497 END_WATCH(BOUNDARIES_EXTERNAL)
00498 #ifdef DEBUG_AMR
00499 CheckLevel(u,Time,Level,t,"AMRSolver<dim>.SetBndry::after SetBndry");
00500 #endif
00501 }
00502
00503 public:
00504 inline const int& Bufferwidth() const { return BufferWidth; }
00505 inline int Bufferwidth() { return BufferWidth; }
00506
00507 protected:
00508 virtual void RegridLevel(const int Time, const int BaseLevel,
00509 BBoxList& bblist, bool TakeListOnFinestLevel) {
00510
00511 #ifdef DEBUG_PRINT_AMRSOLVER
00512 ( comm_service::log() << " *** Regridding at Level: " << BaseLevel
00513 << " t: " << Time << " Max: " << MaxLevel(base::GH()) << " fine: "
00514 << FineLevel(base::GH()) << "\n" ).flush();
00515 #endif
00516
00517 int l;
00518 int flev = FineLevel(base::GH());
00519 if (flev == 0) l = 0;
00520 else if (flev >= MaxLevel(base::GH())) l = MaxLevel(base::GH())-1;
00521 else l = flev;
00522 if (TakeListOnFinestLevel) {
00523 Refine(base::GH(),bblist,l);
00524 l--;
00525 }
00526
00527 for (;l>=BaseLevel;l--) {
00528
00529
00530
00531 START_WATCH
00532
00533 #ifdef DEBUG_AMR
00534 CheckLevel(base::U(),Time,l,t[l],"AMRSolver.RegridLevel::before Flagging");
00535 #endif
00536 base::Flagging_().SetFlags(Time, l, t[l], dt[l]);
00537
00538
00539 forall (base::Flagging_().Flags(),Time,l,c)
00540 for (BBox *_b = bblist.first();_b;_b=bblist.next())
00541 if (!_b->empty()) {
00542 BBox wb = base::Flagging_().Flags().interiorbbox(Time,l,c);
00543 if(wb.stepsize(0)>_b->stepsize(0))
00544 _b->coarsen(wb.stepsize(0)/_b->stepsize(0));
00545 base::Flagging_().Flags()(Time,l,c).plus(200,*_b);
00546 }
00547 end_forall
00548
00549
00550
00551
00552
00553
00554 START_WATCH
00555 Sync(base::Flagging_().Flags(), Time, l);
00556 END_WATCH(BOUNDARIES_SYNC)
00557 END_WATCH_FLAGGING
00558
00559 if (PlotFlags && base::LastOutputTime()==Time) {
00560 char flagname[DAGHBktGFNameWidth+16];
00561 std::sprintf(flagname,"flags");
00562 forall (base::Work(),Time,l,c)
00563 flag_data_type& flag = base::Flagging_().Flags()(Time,l,c);
00564 vflag_grid_data_type vflag_help(flag.bbox(),(VFlagType*) flag.data());
00565 equals_from(base::Work()(Time,l,c), vflag_help, 0);
00566 VFlagType* databuf;
00567 vflag_help.deallocate(databuf);
00568 end_forall
00569 bool flushfile=(l==BaseLevel ? true : false);
00570 base::FileOutput_().WritePlain(base::Work(), flagname, Time, l, t[l], BBox::_empty_bbox, false, flushfile);
00571 }
00572
00573
00574
00575
00576 #ifdef DEBUG_PRINT_INTEGRATOR
00577 ( comm_service::log() << " *** Clustering at Level: " << l << "\n" ).flush();
00578 #endif
00579 START_WATCH
00580 BBox localbbox;
00581 forall(base::Flagging_().Flags(),Time,l,c)
00582 localbbox += base::Flagging_().Flags().boundingbbox(Time,l,c);
00583 end_forall
00584 BBoxList bblexclude, bblcut;
00585 for (BBox *bbi=base::GH().wholebaselist()->first();bbi;
00586 bbi=base::GH().wholebaselist()->next())
00587 bblexclude.add(refine(*bbi,base::GH().refinedby(l)));
00588 bblexclude *= localbbox;
00589 base::GH().glb_bboxlist(bblcut,l,DAGHCellCentered,localbbox);
00590 bblexclude -= bblcut;
00591
00592
00593 DAGHCluster(base::Flagging_().Flags(), Time, l,
00594 BlockWidth, OverlapWidth, BufferWidth,
00595 MinEfficiency, (FlagType)GoodPoint,
00596 bblist, bblist, bblexclude);
00597 END_WATCH_CLUSTERING
00598
00599 if (bblist.isempty()) {
00600 BBox empty_box;
00601 bblist.add(empty_box);
00602 }
00603
00604 #ifdef DEBUG_PRINT_RG_REFINE
00605 ( comm_service::log() << "\n************* New refine list at level "
00606 << l << " *************\n" << bblist
00607 << "\n************* ***************** *************\n"
00608 ).flush();
00609 #endif
00610
00611 Refine(base::GH(),bblist,l);
00612 }
00613 }
00614
00615 virtual void RecomposeGridHierarchy(const int Time, const int Level,
00616 bool ShadowAllowed, bool DoFixup,
00617 bool RecomposeBaseLev, bool RecomposeHighLev) {
00618
00619 if (DoFixup)
00620 base::Fixup_().SetMaxRecomposeLevel(Level);
00621 register int lev;
00622 for (lev=Level; lev<=MaxLevel(base::GH()); lev++) {
00623 base::U().GF_DeleteStorage(1,lev);
00624 if (ShadowAllowed) base::Ush().GF_DeleteStorage(1,lev);
00625 }
00626 #ifdef DEBUG_AMR
00627 if (base::Integrator_().Abort())
00628 for (register int l=0; l<=FineLevel(base::GH()); l++) {
00629 int time = CurrentTime(base::GH(),l);
00630 int tstep = TimeStep(base::U(),l);
00631 CheckLevel(base::U(),time,l,t[l],"AMRSolver.AdvanceLevel::before Recompose");
00632 if (l<Level)
00633 CheckLevel(base::U(),time+tstep,l,t[l],
00634 "AMRSolver.AdvanceLevel::before Recompose for Time+TStep");
00635 if (ShadowAllowed && l<MaxLevel(base::GH())) {
00636 int shtstep = TimeStep(base::Ush(),l);
00637 CheckLevel(base::Ush(),time,l,t[l],
00638 "AMRSolver.AdvanceLevel::before Recompose for Shadow");
00639 if (l<Level)
00640 CheckLevel(base::Ush(),time+shtstep,l,t[l],
00641 "AMRSolver.AdvanceLevel::before Recompose for Time+TStep for Shadow");
00642 }
00643 }
00644 #endif
00645 RecomposeHierarchy(base::GH(), (Level==0 && RecomposeBaseLev) ||
00646 (Level>0 && RecomposeHighLev) ? DAGHTrue : DAGHFalse);
00647 #ifdef DEBUG_AMR
00648 if (base::Integrator_().Abort())
00649 for (register int l=0; l<=FineLevel(base::GH()); l++) {
00650 int time = CurrentTime(base::GH(),l);
00651 int tstep = TimeStep(base::U(),l);
00652 CheckLevel(base::U(),time,l,t[l],"AMRSolver.AdvanceLevel::after Recompose");
00653 if (l<Level)
00654 CheckLevel(base::U(),time+tstep,l,t[l],
00655 "AMRSolver.AdvanceLevel::after Recompose for Time+TStep");
00656 if (ShadowAllowed && l<MaxLevel(base::GH())) {
00657 int shtstep = TimeStep(base::Ush(),l);
00658 CheckLevel(base::Ush(),time,l,t[l],
00659 "AMRSolver.AdvanceLevel::after Recompose for Shadow");
00660 if (l<Level)
00661 CheckLevel(base::Ush(),time+shtstep,l,t[l],
00662 "AMRSolver.AdvanceLevel::after Recompose for Time+TStep for Shadow");
00663 }
00664 }
00665 #endif
00666
00667 for (lev=Level; lev<=MaxLevel(base::GH()); lev++) {
00668 base::U().GF_CreateStorage(1,lev);
00669 if (ShadowAllowed && lev<MaxLevel(base::GH()))
00670 base::Ush().GF_CreateStorage(1,lev);
00671 }
00672 }
00673
00674 virtual void BeforeLevelStep(const int Level) {}
00675 virtual void AfterLevelStep(const int Level) {}
00676
00677 virtual void AdvanceLevel(const int Level, int RegridEvery, bool RegridDone,
00678 bool ShadowAllowed, bool DoFixup,
00679 bool RecomposeBaseLev, bool RecomposeHighLev) {
00680
00681
00682
00683 int NoIterations = 0;
00684 if (Level==0 || NoTimeRefine==1) {
00685 NoIterations = 1;
00686 dt[Level] = DeltaT(base::GH(), 0);
00687 }
00688 else {
00689 NoIterations = RefineFactor(base::GH(),Level-1);
00690 dt[Level] = dt[Level-1] / RefineFactor(base::GH(),Level-1);
00691 }
00692
00693 for (int l=Level+1; l<=MaxLevel(base::GH()); l++)
00694 if (NoTimeRefine==1) dt[l]=dt[l-1];
00695 else dt[l] = dt[l-1] / RefineFactor(base::GH(),l-1);
00696
00697 int Time;
00698 int TStep = TimeStep(base::U(),Level);
00699 double cfl;
00700 cfl_new[Level] = 0.0;
00701
00702 for (int k=0; k<NoIterations; k++) {
00703 Time = CurrentTime(base::GH(),Level);
00704
00705 #ifdef DEBUG_PRINT_AMRSOLVER
00706 ( comm_service::log() << "AdvanceLevel(): Level" << Level
00707 << " Time = " << Time << " Iteration = " << k << "\n" ).flush();
00708 #endif
00709
00710 SetPhysicalTime(base::U(),Time+TStep,Level,t[Level]+dt[Level]);
00711 if (ShadowAllowed)
00712 SetPhysicalTime(base::Ush(),Time+TimeStep(base::Ush(),Level),Level,
00713 t[Level]+dt[Level]*base::Ush().factor());
00714
00715 bool DoRegrid = false;
00716
00717 if (NoTimeRefine==1) {
00718 if (Level==0)
00719 for (int l=0; l<=FineLevel(base::GH()); l++)
00720 SetBndry(base::U(),Time,l,t[l]);
00721
00722 DoRegrid = (RegridEvery!= AMRSolverRegridFalse ?
00723 Level==0 && StepsTaken(base::GH(),Level)%RegridEvery==0 &&
00724 t[Level]>0.0 : false);
00725 }
00726 else {
00727 if (Level==0 || !RegridDone)
00728 SetBndry(base::U(),Time,Level,t[Level]);
00729
00730 if (RegridEvery>AMRSolverRegridFalse)
00731 DoRegrid = Level<MaxLevel(base::GH()) && k%RegridEvery==0 && !RegridDone && t[Level]>0.0;
00732 else if (RegridEvery<AMRSolverRegridFalse)
00733 DoRegrid = Level==0 && StepsTaken(base::GH(),Level)%RegridEvery==0 && t[Level]>0.0;
00734 else
00735 DoRegrid = false;
00736 }
00737
00738
00739
00740
00741 if (DoRegrid) {
00742
00743
00744
00745 BBoxList bblist;
00746 RegridLevel(Time, Level, bblist, false);
00747
00748
00749
00750
00751 START_WATCH
00752 RecomposeGridHierarchy(Time,Level,ShadowAllowed,DoFixup,
00753 RecomposeBaseLev,RecomposeHighLev);
00754 END_WATCH_RECOMPOSING_WHOLE
00755 RegridDone = true;
00756 }
00757
00758 BeforeLevelStep(Level);
00759
00760
00761
00762
00763
00764
00765
00766
00767
00768 cfl = IntegrateLevel(base::U(), Time, Level, t[Level], dt[Level],
00769 DoFixup, (Level>0 ? t[Level-1] : 0.0), 0);
00770 cfl_new[Level] = cfl_new[Level] > cfl ? cfl_new[Level] : cfl;
00771
00772
00773
00774
00775
00776
00777
00778
00779 if (RegridEvery != AMRSolverRegridFalse ?
00780 Level<MaxLevel(base::GH()) && ShadowAllowed &&
00781 (Level>0 ? (k+1)%RegridEvery==0 : true) : false) {
00782
00783 double dt_sh = dt[Level] * base::Ush().factor();
00784 SetPhysicalTime(base::Ush(),Time,Level,t[Level]);
00785
00786
00787
00788
00789 forall (base::U(),Time,Level,cm)
00790 vec_grid_data_type& u_main = base::U()(Time, Level, cm);
00791 forall (base::Ush(),Time,Level,cs)
00792 vec_grid_data_type& u_shadow = base::Ush()(Time, Level, cs);
00793 BBox bb = u_shadow.bbox()*u_main.bbox();
00794 if (!bb.empty())
00795 base::LevelTransfer_().Restrict(u_main,Level,u_shadow,Level,bb);
00796 end_forall
00797 end_forall
00798
00799 if (NoTimeRefine!=1)
00800 SetBndry(base::Ush(),Time,Level,t[Level]);
00801 IntegrateLevel(base::Ush(), Time, Level, t[Level], dt_sh, false, 0.0, 2);
00802
00803 }
00804
00805
00806
00807
00808
00809 #ifdef DEBUG_PRINT_FIXUP
00810 VectorType mass_old;
00811 int d;
00812 DCoords dx = base::GH().worldStep(base::U().GF_StepSize(Level));
00813 if (Level == 0) {
00814 mass_old = Sum(base::U(),Time+TStep,Level);
00815 for (d=0; d<dim; d++) mass_old *= dx(d);
00816 }
00817 #endif
00818
00819 if (Level < MaxLevel(base::GH())) {
00820
00821
00822
00823
00824 if (NoTimeRefine!=1 && AdaptBndTimeInterpolate!=0)
00825 SetBndry(base::U(),Time+TStep,Level,t[Level]+dt[Level]);
00826
00827
00828 AdvanceLevel(Level+1, RegridEvery, RegridDone,
00829 ShadowAllowed, DoFixup, RecomposeBaseLev,
00830 RecomposeHighLev);
00831
00832 #ifdef DEBUG_AMR
00833 if (AdaptBndTimeInterpolate!=0)
00834 CheckLevel(base::U(),Time+TStep,Level,t[Level],
00835 "AMRSolver.AdvanceLevel::before Restrict and Fixup");
00836 #endif
00837
00838
00839
00840 if (DoFixup) {
00841 START_WATCH
00842 base::Fixup_().Correction(Time+TStep, Time, Level, t[Level], dt[Level]);
00843 END_WATCH(FIXUP_CORRECTION)
00844 }
00845 int TimeF = CurrentTime(base::GH(),Level+1);
00846 Restrict(base::U(), TimeF, Level+1, Time+TStep, Level);
00847
00848 #ifdef DEBUG_AMR
00849 if (AdaptBndTimeInterpolate!=0)
00850 CheckLevel(base::U(),Time+TStep,Level,t[Level],
00851 "AMRSolver.AdvanceLevel::after Restrict and Fixup");
00852 #endif
00853 }
00854 #ifdef DEBUG_PROFILE
00855 ( std::cout << "\n" ).flush();
00856 #endif
00857
00858 #ifdef DEBUG_PRINT_FIXUP
00859 VectorType mass_new;
00860 if (Level == 0) {
00861 ( comm_service::log()
00862 << "\n*********** Difference in state variables at level "
00863 << Level << " after fixup ********\n" ).flush();
00864
00865
00866 mass_new = Sum(base::U(),Time+TStep,Level);
00867 for (d=0; d<dim; d++) mass_new *= dx(d);
00868 ( comm_service::log() << mass_old << "-"
00869 << mass_new << "=" << mass_old - mass_new ).flush();
00870 ( comm_service::log()
00871 << "\n**************************************************************"
00872 << "***********\n" ).flush();
00873 }
00874 #endif
00875
00876
00877
00878
00879 CycleTimeLevels(base::U(),Level);
00880 if (ShadowAllowed)
00881 CycleTimeLevels(base::Ush(),Level);
00882
00883
00884
00885
00886
00887
00888
00889 t[Level] += dt[Level];
00890 IncrCurrentTime(base::GH(),Level);
00891 if (NoTimeRefine==1 && Level>0)
00892 for (int kk=1; kk<RefinedBy(base::GH(),Level); kk++)
00893 IncrCurrentTime(base::GH(),Level);
00894 Time = CurrentTime(base::GH(),Level);
00895 SetPhysicalTime(base::U(),Time,Level,t[Level]);
00896 if (ShadowAllowed)
00897 SetPhysicalTime(base::Ush(),Time,Level,t[Level]+
00898 dt[Level]*(base::Ush().factor()-1));
00899 RegridDone = false;
00900
00901 AfterLevelStep(Level);
00902 }
00903 }
00904
00905 void CheckLevel(vec_grid_fct_type &u, const int& time, const int& level,
00906 const double& t, const char *text) {
00907 if (base::Integrator_().Abort()) {
00908 char wtext[1024];
00909 std::sprintf(wtext,"\non Level=%d at Time=%d: %s\n",level,time,text);
00910 forall (u,time,level,c)
00911 vec_grid_data_type& uw = u(time,level,c);
00912 base::Integrator_().CheckGrid(uw, level, uw.bbox(), t, wtext);
00913 end_forall
00914 }
00915 }
00916
00917 BBox BoundaryBBox(const BBox& whole, const int s) {
00918 int d = s/2;
00919 #ifdef DEBUG_PRINT
00920 assert (d>=0 && d<dim);
00921 #endif
00922 BBox bb(whole);
00923 if (s%2 == 0)
00924 bb.upper(d) = bb.lower(d)+(base::NGhosts()-1)*bb.stepsize(d);
00925 else
00926 bb.lower(d) = bb.upper(d)-(base::NGhosts()-1)*bb.stepsize(d);
00927 return bb;
00928 }
00929
00930 int FixupPar;
00931 int RegridEvery;
00932 double MinEfficiency;
00933 int BlockWidth, OverlapWidth, BufferWidth, NestingBuffer;
00934 int NoTimeRefine, AdaptBndTimeInterpolate;
00935 int PlotFlags;
00936
00937 double *t, *dt;
00938 double *cfl_new;
00939 };
00940
00941
00942 #endif