Blockstructured Adaptive Mesh Refinement in object-oriented C++
Blockstructured Adaptive Mesh Refinement in object-oriented C++
00001 #ifndef AMROC_AMRSOLVER_H 00002 #define AMROC_AMRSOLVER_H 00003 00011 #include "Solver.h" 00012 #include "SolverControl.h" 00013 #include "AMRSolverInterface.h" 00014 00015 00016 #ifndef AMRSolverName 00017 #define AMRSolver(dim) name2(AMRSolver,dim) 00018 #define AMRSolverName 00019 #endif 00020 00021 #define DAGHCluster(dim) name2(name2(DAGHCluster,dim),d) 00022 00023 #define AMRSolverRegridFalse (0) 00024 #define VERY_LARGE_CFL (1.5) 00025 #define MAXOVERLAPWIDTH (100) 00026 00033 //****************************************************************************** 00034 // flag_cells_by_error(), flag_cells_by_difference() and the fixup make use of 00035 // u(Time+TStep). Ensure that the next time step is only used when all data, 00036 // especially the adaptive boundary conditions (-> RecomposeHierarchy() ) have 00037 // been set. 00038 // If problems occur, the use of a working function might be unavoidable. 00039 //****************************************************************************** 00040 00041 template <class VectorType, class FixupType, class FlagType> 00042 class AMRSolver(DIM) : 00043 public Solver(DIM)<VectorType,FixupType,FlagType>, 00044 public AMRSolverInterface(DIM)<VectorType> { 00045 00046 typedef Integrator(DIM)<VectorType> integrator_type; 00047 typedef InitialCondition(DIM)<VectorType> initial_condition_type; 00048 typedef BoundaryConditions(DIM)<VectorType> boundary_conditions_type; 00049 00050 typedef typename VectorType::InternalDataType DataType; 00051 typedef GridFunction(DIM)<VectorType> vec_grid_fct_type; 00052 typedef GridData(DIM)<VectorType> vec_grid_data_type; 00053 00054 typedef Vector(1)<FlagType> VFlagType; 00055 typedef GridData(DIM)<VFlagType> vflag_grid_data_type; 00056 typedef GridData(DIM)<FlagType> flag_data_type; 00057 00058 public: 00059 AMRSolver(DIM)(integrator_type& integ, 00060 initial_condition_type& init, 00061 boundary_conditions_type& bc) : 00062 Solver(DIM)<VectorType, FixupType, FlagType>(integ, init, bc), 00063 FixupPar(1), RegridEvery(1), MinEfficiency(0.7), 00064 BlockWidth(1), OverlapWidth(0), BufferWidth(2), NestingBuffer(1), 00065 PlotFlags(0) { 00066 00067 t = (double *)0; dt = (double *)0; cfl_new = (double *)0; 00068 BufferWidth = NGhosts(); 00069 AllocError::SetTexts("AMRSolver(DIM)","Constructor"); 00070 } 00071 00072 ~AMRSolver(DIM)() { 00073 if (t) delete [] t; 00074 if (dt) delete [] dt; 00075 if (cfl_new) delete [] cfl_new; 00076 } 00077 00078 virtual void register_at(ControlDevice& Ctrl, const string& prefix) { 00079 LocCtrl = Ctrl.getSubDevice(prefix+"AMRSolver"); 00080 RegisterAt(LocCtrl,"DoFixup",FixupPar); 00081 RegisterAt(LocCtrl,"RegridEvery",RegridEvery); 00082 RegisterAt(LocCtrl,"MinEfficiency",MinEfficiency); 00083 RegisterAt(LocCtrl,"BufferWidth",BufferWidth); 00084 RegisterAt(LocCtrl,"BlockWidth",BlockWidth); 00085 RegisterAt(LocCtrl,"OverlapWidth",OverlapWidth); 00086 RegisterAt(LocCtrl,"NestingBuffer",NestingBuffer); 00087 RegisterAt(LocCtrl,"PlotFlags",PlotFlags); 00088 00089 Integrator_().register_at(LocCtrl, prefix); 00090 InitialCondition_().register_at(LocCtrl, prefix); 00091 BoundaryConditions_().register_at(LocCtrl, prefix); 00092 if (_Flagging) 00093 Flagging_().register_at(LocCtrl, prefix); 00094 if (_Fixup) 00095 Fixup_().register_at(LocCtrl, prefix); 00096 } 00097 virtual void register_at(ControlDevice& Ctrl) { 00098 register_at(Ctrl, ""); 00099 } 00100 00101 virtual void SetupData() { 00102 int me = MY_PROC(_Hierarchy); 00103 int help_set = NGhosts(); 00104 f_init_opcommon(help_set); 00105 SetProlongFunction(U(), (void *) &f_prolong_amr); 00106 SetRestrictFunction(U(), (void *) &f_restrict_amr); 00107 if (ErrorEstimation()) { 00108 SetProlongFunction(Ush(), (void *) &f_prolong_amr); 00109 SetRestrictFunction(Ush(), (void *) &f_restrict_amr); 00110 } 00111 00112 t = new double[MaxLevel(GH())+1]; 00113 dt = new double[MaxLevel(GH())+1]; 00114 cfl_new = new double[MaxLevel(GH())+1]; 00115 00116 Integrator_().SetupData(); 00117 if (_Flagging && MaxLevel(GH())>=1) { 00118 Flagging_().SetBufferwidth(Min(Max(BufferWidth,OverlapWidth-1), 00119 MAXOVERLAPWIDTH)); 00120 Flagging_().SetupData(); 00121 } 00122 else 00123 RegridEvery = 0; 00124 00125 if (_Fixup && FixupPar>0) 00126 Fixup_().SetupData(); 00127 else 00128 FixupPar = 0; 00129 00130 if (FixupPar && NestingBuffer<1 && RegridEvery!=AMRSolverRegridFalse) { 00131 if (me == VizServer) 00132 cout << " Fixup requires a NestingBuffer>0 ! Using Default." << endl; 00133 NestingBuffer = 1; 00134 } 00135 GH().DAGH_SetNestingBuffer(NestingBuffer); 00136 } 00137 00138 virtual void Initialize(const double dt_start) { 00139 START_WATCH 00140 00141 U().GF_DeleteStorage(1); 00142 if (ErrorEstimation()) 00143 Ush().GF_DeleteStorage(1); 00144 00145 //--------------------------------------------------------------- 00146 // Calculate first grid distribution 00147 //--------------------------------------------------------------- 00148 int lev = 0; 00149 if (RegridEvery != AMRSolverRegridFalse) { 00150 while (lev<=FineLevel(GH()) && lev<MaxLevel(GH())) { 00151 //------------------------------------------------------------------ 00152 // Set up initial data 00153 //------------------------------------------------------------------ 00154 #ifdef DEBUG_PRINT_AMRSOLVER 00155 ( comm_service::log() << " *** Setting flags at Level: " << lev << "\n" ).flush(); 00156 #endif 00157 //------------------------------------------------------------------ 00158 // Initializing all levels ensures correct prolongation 00159 //------------------------------------------------------------------ 00160 for (register int l=0; l<=FineLevel(GH()); l++) { 00161 InitialCondition_().Set(U(), l); 00162 if (ErrorEstimation()) 00163 InitialCondition_().Set(Ush(), l); 00164 } 00165 BBoxList bblist; 00166 RegridLevel(0, 0, bblist, false); 00167 00168 END_WATCH_INITIALIZATION 00169 START_WATCH 00170 RecomposeHierarchy(GH()); 00171 END_WATCH_RECOMPOSING_WHOLE 00172 START_WATCH 00173 lev++; 00174 } 00175 00176 if (FineLevel(GH())>0) { 00177 for (register int l=0; l<=FineLevel(GH()); l++) 00178 InitialCondition_().Set(U(), l); 00179 BBoxList bblist; 00180 forall (U(),0,FineLevel(GH()),c) 00181 bblist.add(coarsen(U().interiorbbox(0,FineLevel(GH()),c), 00182 RefineFactor(GH(),FineLevel(GH())-1))); 00183 end_forall 00184 RegridLevel(0, 0, bblist, true); 00185 END_WATCH_INITIALIZATION 00186 START_WATCH 00187 RecomposeHierarchy(GH()); 00188 END_WATCH_RECOMPOSING_WHOLE 00189 START_WATCH 00190 } 00191 } 00192 00193 U().GF_CreateStorage(1); 00194 if (ErrorEstimation()) 00195 Ush().GF_CreateStorage(1); 00196 00197 for (lev=0; lev<=FineLevel(GH()); lev++) { 00198 #ifdef DEBUG_PRINT_AMRSOLVER 00199 ( comm_service::log() << " *** Initializing Level: " << lev << " *** \n" ).flush(); 00200 #endif 00201 InitialCondition_().Set(U(), lev); 00202 if (ErrorEstimation()) 00203 InitialCondition_().Set(Ush(), lev); 00204 t[lev] = 0.0; 00205 } 00206 00207 SetTimeSpecs(GH(), 0.0, dt_start, 2); 00208 END_WATCH_INITIALIZATION 00209 } 00210 00211 virtual void Restart(const char* CheckpointFile) { 00212 START_WATCH 00213 RecomposeHierarchy(GH(), CheckpointFile); 00214 double NextTime; 00215 for (register int l=0; l<=FineLevel(GH()); l++) { 00216 GH().DAGH_GetTimeSpecs(t[l], NextTime); 00217 int Time = CurrentTime(GH(),l); 00218 SetPhysicalTime(U(),Time,l,t[l]); 00219 if (ErrorEstimation()) SetPhysicalTime(Ush(),Time,l,t[l]); 00220 } 00221 END_WATCH_INITIALIZATION 00222 } 00223 00224 virtual double Tick(int VariableTimeStepping, const double dtv[], const double cflv[], 00225 int& Rejections, int& RedistributeEvery) { 00226 00227 double told = t[0]; 00228 double cfl_max; 00229 Rejections = 0; 00230 do { 00231 strstream* RestartStorage = (strstream *)0; 00232 char* RestartPt = (char *)0; 00233 if (VariableTimeStepping==AutomaticTimeStepsRestart) { 00234 AllocError::SetTexts("AMRSolver(DIM)","Memory for restart"); 00235 #if defined(__xlC__) 00236 RestartStorage = new strstream(); 00237 #else 00238 int sizept = GH().DAGH_Checkpoint_StrStream_Memory(); 00239 RestartPt = new char[sizept]; 00240 RestartStorage = new strstream(RestartPt,sizept); 00241 #endif 00242 GH().DAGH_Checkpoint(*RestartStorage); 00243 } 00244 00245 #ifdef DEBUG_PRINT_AMRSOLVER 00246 double Current_Time, Next_Time; 00247 GH().DAGH_GetTimeSpecs(Current_Time, Next_Time); 00248 ( comm_service::log() << " --- Iteration: " << StepsTaken(GH(),0)+1 << " ---" 00249 << " t = " << Current_Time << " dt = " 00250 << DeltaT(GH(), 0) << "\n" ).flush(); 00251 #endif 00252 #ifdef DEBUG_PROFILE 00253 ( cout << "\n" ).flush(); 00254 #endif 00255 00256 dt[0] = DeltaT(GH(), 0); 00257 int l=1; 00258 for (; l<=FineLevel(GH()); l++) 00259 dt[l] = dt[l-1] / RefineFactor(GH(),l-1); 00260 // Use this for simplified AMR 00261 // dt[l] = dt[l-1]; 00262 00263 AdvanceLevel(0, RegridEvery, 00264 !(RegridEvery != AMRSolverRegridFalse ? 00265 StepsTaken(GH(),0)%RegridEvery==0 && StepsTaken(GH(),0)>0 : false), 00266 ErrorEstimation() && (RegridEvery != AMRSolverRegridFalse ? 00267 (StepsTaken(GH(),0)+1)%RegridEvery==0 : false), FixupPar > 0, 00268 !(RedistributeEvery>0 && StepsTaken(GH(),0)%RedistributeEvery!=0), 00269 RedistributeEvery<=0); 00270 00271 cfl_max = 0.0; 00272 for (l=0; l<=FineLevel(GH()); l++) 00273 cfl_max = cfl_max > cfl_new[l] ? cfl_max : cfl_new[l]; 00274 00275 #ifdef DAGH_NO_MPI 00276 #else 00277 if (comm_service::dce() && comm_service::proc_num() > 1) { 00278 00279 int num = comm_service::proc_num(); 00280 00281 int R; 00282 int size = sizeof(double); 00283 void *values = (void *) new char[size*num]; 00284 00285 R = MPI_Allgather(&cfl_max, size, MPI_BYTE, values, size, MPI_BYTE, 00286 comm_service::comm(U().GF_Id())); 00287 if ( MPI_SUCCESS != R ) 00288 comm_service::error_die( "AMRSolver(DIM)::cfl_max", "MPI_Allgather", R ); 00289 00290 for (int i=0;i<num;i++) { 00291 double tmax = *((double *)values + i); 00292 cfl_max = cfl_max > tmax ? cfl_max : tmax; 00293 } 00294 if (values) delete [] (char *)values; 00295 } 00296 #endif 00297 00298 #ifdef DEBUG_PRINT_AMRSOLVER 00299 ( comm_service::log() << " Levels = " << FineLevel(GH())+1 00300 << " max. CFL = " << cfl_max << "\n" ).flush(); 00301 #endif 00302 #ifdef DEBUG_PROFILE 00303 ( cout << "\n" ).flush(); 00304 #endif 00305 00306 if (cfl_max > cflv[0]) 00307 if (VariableTimeStepping==AutomaticTimeStepsRestart) { 00308 if (RestartStorage!=(strstream *)0) { 00309 GH().DAGH_RecomposeHierarchy(*RestartStorage); 00310 for (int lev=0; lev<=FineLevel(GH()); lev++) 00311 t[lev] = told; 00312 Rejections++; 00313 } 00314 else 00315 cerr << " Restart failure!" << endl << flush ; 00316 if (cfl_max > VERY_LARGE_CFL) 00317 cerr << " max. CFL>" << (double)VERY_LARGE_CFL << " " << flush; 00318 } 00319 else if (cfl_max > 1.0) 00320 cerr << " max. CFL>1.0 ! " << flush; 00321 00322 if (VariableTimeStepping!=FixedTimeSteps) 00323 if (cfl_max == 0.0) 00324 SetTimeSpecs(GH(), t[0], t[0]+dtv[1], 2); 00325 else { 00326 double dt_predict = (DeltaT(GH(), 0) / cfl_max) * cflv[1]; 00327 SetTimeSpecs(GH(), t[0], t[0] + 00328 (dtv[1] < dt_predict ? dtv[1] : dt_predict), 2); 00329 } 00330 else 00331 SetTimeSpecs(GH(), t[0], t[0]+dtv[0], 2); 00332 00333 if (RestartStorage != (strstream *)0) 00334 delete RestartStorage; 00335 if (RestartPt != (char *)0) 00336 delete [] RestartPt; 00337 00338 } while (t[0] <= told); 00339 00340 return cfl_max; 00341 } 00342 00343 virtual double IntegrateLevel(vec_grid_fct_type& u, const int Time, const int Level, 00344 double t, double dt, bool DoFixup, double tc, 00345 const int which) { 00346 00347 //--------------------------------------------------------------- 00348 // Take one step on u. u(Time, Level) remains unchanged. 00349 //--------------------------------------------------------------- 00350 char errtext[256], whichtext[16]; 00351 if (which<=0 || which>2) sprintf(whichtext,"Main"); 00352 else if (which==1) sprintf(whichtext,"Estimate"); 00353 else sprintf(whichtext,"Shadow"); 00354 sprintf(errtext,"Flux allocation for %s",whichtext); 00355 AllocError::SetTexts("AMRSolver(DIM)",errtext); 00356 00357 int TStep = TimeStep(u,Level); 00358 double cfl = 0, cfl_grid; 00359 DCoords dx(DIM); 00360 for (int d=0; d<DIM; d++) 00361 dx(d) = DeltaX(u,d,Level); 00362 00363 int mpass = 0; 00364 if (Integrator_().MaxIntegratorPasses() > 0) mpass = 1; 00365 do { 00366 00367 forall (u,Time,Level,c) 00368 vec_grid_data_type& u_old = u(Time, Level, c); 00369 vec_grid_data_type& u_new = u(Time+TStep, Level, c); 00370 00371 temporary_fluxes(f, u_old) 00372 #ifdef DEBUG_PRINT_INTEGRATOR 00373 if (mpass<=1) ( comm_service::log() << " " << which << ": " ).flush(); 00374 #endif 00375 START_WATCH 00376 cfl_grid = Integrator_().CalculateGrid(u_new, u_old, f, t, dt, dx, mpass); 00377 END_WATCH_INTEGRATION(which) 00378 #ifdef DEBUG_PROFILE 00379 if (mpass<=1) ( cout << "." ).flush(); 00380 #endif 00381 00382 cfl = (cfl > cfl_grid ? cfl : cfl_grid); 00383 00384 if (Level>0 && DoFixup) { 00385 // add_first_order_fluxes() solves a Riemann-problem on the basis of 00386 // u(CurrentTime(GH,Level-1),Level-1) at fine-coarse interfaces. 00387 // Do not change u(CurrentTime(GH,Level-1),Level-1)! 00388 START_WATCH 00389 Fixup_().AddFluxes(Time, Level, c, f, tc, t, dt, dx, mpass); 00390 END_WATCH_FIXUP_WHOLE 00391 } 00392 00393 if (Level<FineLevel(GH()) && DoFixup) { 00394 START_WATCH 00395 Fixup_().SaveFluxes(Time, Level, c, f , dt, dx, mpass); 00396 END_WATCH_FIXUP_WHOLE 00397 } 00398 00399 end_temporary_fluxes(f) 00400 end_forall 00401 00402 // Set boundary values again at t=Time. As the partially updated time 00403 // step is stored at Time+TStep the two time steps have to swaped. 00404 if (mpass > 0) { 00405 SwapTimeLevels(u,Level,Time,Time+TStep); 00406 START_WATCH 00407 AdaptiveBoundaryUpdate(u,Time,Level); 00408 END_WATCH_BOUNDARIES_INTERPOLATION 00409 START_WATCH 00410 Sync(u,Time,Level); 00411 END_WATCH_BOUNDARIES_SYNC 00412 START_WATCH 00413 ExternalBoundaryUpdate(u,Time,Level); 00414 END_WATCH_BOUNDARIES_WHOLE 00415 SwapTimeLevels(u,Level,Time,Time+TStep); 00416 } 00417 00418 mpass++; 00419 } while (mpass <= Integrator_().MaxIntegratorPasses()); 00420 00421 return cfl; 00422 } 00423 00424 protected: 00425 void RegridLevel(const int Time, const int BaseLevel, 00426 BBoxList& bblist, bool TakeListOnFinestLevel) { 00427 00428 #ifdef DEBUG_PRINT_AMRSOLVER 00429 ( comm_service::log() << " *** Regridding at Level: " << BaseLevel 00430 << " t: " << Time << " Max: " << MaxLevel(GH()) << " fine: " 00431 << FineLevel(GH()) << "\n" ).flush(); 00432 #endif 00433 00434 int l; 00435 int flev = FineLevel(GH()); 00436 if (flev == 0) l = 0; 00437 else if (flev >= MaxLevel(GH())) l = MaxLevel(GH())-1; 00438 else l = flev; 00439 if (TakeListOnFinestLevel) { 00440 Refine(GH(),bblist,l); 00441 l--; 00442 } 00443 00444 char flagname[DAGHBktGFNameWidth+16]; 00445 if (PlotFlags && LastOutputTime()==Time) 00446 sprintf(flagname,"flag_%d",Time); 00447 int me = MY_PROC(GH); 00448 00449 for (;l>=BaseLevel;l--) { 00450 //--------------------------------------------------------------- 00451 // Flag cells for refinement by error estimation 00452 //--------------------------------------------------------------- 00453 START_WATCH 00454 00455 #ifdef DEBUG_AMR 00456 CheckLevel(U(),Time,l,"AMRSolver(DIM).RegridLevel::before Flagging"); 00457 #endif 00458 Flagging_().FlagLevel(Time, l, t[l], dt[l]); 00459 00460 // Ensure nesting of higher level patches. Then syncronize the flags. 00461 forall (Flagging_().Flags(),Time,l,c) 00462 for (BBox *_b = bblist.first();_b;_b=bblist.next()) 00463 if (!_b->empty()) { 00464 BBox wb = Flagging_().Flags().interiorbbox(Time,l,c); 00465 if(wb.stepsize(0)>_b->stepsize(0)) 00466 _b->coarsen(wb.stepsize(0)/_b->stepsize(0)); 00467 Flagging_().Flags()(Time,l,c).plus(200,*_b); 00468 } 00469 end_forall 00470 00471 // The flags have to be synced to create the bufferzone around flagged cells 00472 // correctly. Without syncing no bufferzone into neighbouring patches on other 00473 // processors would be created. 00474 // The radius for the syncronization stencil is set to 00475 // Min(Max(BufferWidth,OverlapWidth-1),MAXOVERLAPWIDTH) 00476 START_WATCH 00477 Sync(Flagging_().Flags(), Time, l); 00478 END_WATCH_BOUNDARIES_SYNC 00479 END_WATCH_FLAGGING 00480 00481 if (PlotFlags && LastOutputTime()==Time) { 00482 forall (Work(),Time,l,c) 00483 flag_data_type& flag = Flagging_().Flags()(Time,l,c); 00484 vflag_grid_data_type vflag_help(flag.bbox(),(VFlagType*) flag.data()); 00485 equals_from(Work()(Time,l,c), vflag_help, 0); 00486 VFlagType* databuf; 00487 vflag_help.deallocate(databuf); 00488 end_forall 00489 if (me == VizServer && l==0) 00490 cout << " *** Writing " << flagname << endl; 00491 Write(Work(),Time,l,DAGH_Double,flagname); 00492 } 00493 00494 //--------------------------------------------------------------- 00495 // Cluster at refine 00496 //--------------------------------------------------------------- 00497 #ifdef DEBUG_PRINT_INTEGRATOR 00498 ( comm_service::log() << " *** Clustering at Level: " << l << "\n" ).flush(); 00499 #endif 00500 START_WATCH 00501 BBoxList bblexclude, bblcut; 00502 for (BBox *bbi=GH().wholebaselist()->first();bbi; 00503 bbi=GH().wholebaselist()->next()) 00504 bblexclude.add(refine(*bbi,GH().refinedby(l))); 00505 GH().glb_bboxlist(bblcut,l); 00506 bblexclude -= bblcut; 00507 00508 // After syncing, the clusterer creates the bufferzone locally. 00509 DAGHCluster(DIM)(Flagging_().Flags(), Time, l, 00510 BlockWidth, OverlapWidth, BufferWidth, 00511 MinEfficiency, (FlagType)GoodPoint, 00512 bblist, bblist, bblexclude); 00513 END_WATCH_CLUSTERING 00514 00515 if (bblist.isempty()) { 00516 BBox empty_box; 00517 bblist.add(empty_box); 00518 } 00519 00520 #ifdef DEBUG_PRINT_RG_REFINE 00521 ( comm_service::log() << "\n************* New refine list at level " 00522 << l << " *************\n" << bblist 00523 << "\n************* ***************** *************\n" 00524 ).flush(); 00525 #endif 00526 00527 Refine(GH(),bblist,l); 00528 } 00529 00530 #ifndef DAGH_NO_MPI 00531 if (PlotFlags && LastOutputTime()==Time) 00532 MPI_Barrier(comm_service::comm()); 00533 #endif 00534 } 00535 00536 virtual void AdvanceLevel(const int Level, int RegridEvery, bool RegridDone, 00537 bool ShadowAllowed, bool DoFixup, 00538 bool RecomposeBaseLev, bool RecomposeHighLev) { 00539 //--------------------------------------------------------------- 00540 // Select number of iterations to run based on current grid level 00541 //--------------------------------------------------------------- 00542 int NoIterations = 0; 00543 if (Level==0) 00544 NoIterations = 1; 00545 else 00546 NoIterations = RefineFactor(GH(),Level-1); 00547 // Use this for simplified AMR 00548 // NoIterations = 1; 00549 00550 int Time; 00551 int TStep = TimeStep(U(),Level); 00552 double cfl; 00553 cfl_new[Level] = 0.0; 00554 DCoords dx(DIM); 00555 for(int d=0; d<DIM; d++) 00556 dx(d) = DeltaX(U(),d,Level); 00557 00558 for (int k=0; k<NoIterations; k++) { 00559 Time = CurrentTime(GH(),Level); 00560 SetPhysicalTime(U(),Time+TStep,Level,t[Level]+dt[Level]); 00561 if (ShadowAllowed) 00562 SetPhysicalTime(Ush(),Time+TimeStep(Ush(),Level),Level, 00563 t[Level]+dt[Level]*Ush().factor()); 00564 00565 bool DoRegrid = (RegridEvery!= AMRSolverRegridFalse ? 00566 Level<MaxLevel(GH()) && k%RegridEvery==0 && 00567 !RegridDone && t[Level]>0.0 : false); 00568 00569 // Set boundary conditions for current time 00570 // 00571 // Interpolation at fine-coarse interfaces is done in the following way: 00572 // k=0: Space prolongation from 00573 // u(CurrentTime(GH,Level-1),Level-1) 00574 // k>0: Space prolongation of 00575 // u(CurrentTime(GH,Level-1)+TimeStep(GH,Level-1), 00576 // Level-1) and (CurrentTime(GH,Level-1,) and 00577 // time interpolation 00578 // Ghostcells are filled in the following order: 00579 // 1. forall (u,Time,Level,c) 00580 // Interpolate completely into all sides abutting a patch of the coarser level. 00581 // 2. forall (u,Time,Level,c) 00582 // Override previous values in ghostcells overlying a patch of the same level. 00583 // 3. forall (u,Time,Level,c) 00584 // Set physical boundary conditions at sides outside of the domain. 00585 // BoundaryUpdate() calls all necessary functions in the right order! 00586 // RecomposeHierarchy() calls BoundaryUpdate() on all available time steps 00587 // and levels. Therefore, we can skip the setting of boundary conditions on fine 00588 // grid levels, if RecomposeHierarchy() has been called on the coarse base level. 00589 if (Level==0 || !RegridDone) { 00590 START_WATCH 00591 AdaptiveBoundaryUpdate(U(),Time,Level); 00592 #ifdef DEBUG_AMR 00593 CheckLevel(U(),Time,Level, 00594 "AMRSolver(DIM).AdvanceLevel::after AdaptiveBoundaryUpdate"); 00595 #endif 00596 END_WATCH_BOUNDARIES_INTERPOLATION 00597 00598 START_WATCH 00599 Sync(U(),Time,Level); 00600 #ifdef DEBUG_AMR 00601 CheckLevel(U(),Time,Level,"AMRSolver(DIM).AdvanceLevel::after Sync"); 00602 #endif 00603 END_WATCH_BOUNDARIES_SYNC 00604 00605 START_WATCH 00606 ExternalBoundaryUpdate(U(),Time,Level); 00607 #ifdef DEBUG_AMR 00608 CheckLevel(U(),Time,Level, 00609 "AMRSolver(DIM).AdvanceLevel::after ExternalBoundaryUpdate"); 00610 #endif 00611 END_WATCH_BOUNDARIES_WHOLE 00612 } 00613 00614 //--------------------------------------------------------------- 00615 // regridding... 00616 //--------------------------------------------------------------- 00617 if (DoRegrid) { 00618 //--------------------------------------------------------------- 00619 // Regrid all levels > Level, Level itself stays fixed 00620 //--------------------------------------------------------------- 00621 BBoxList bblist; 00622 RegridLevel(Time, Level, bblist, false); 00623 00624 //--------------------------------------------------------------- 00625 // Recompose the GridHierarchy 00626 //--------------------------------------------------------------- 00627 int finest = FineLevel(GH()); 00628 00629 START_WATCH 00630 // Time interpolation needs two time steps only on the coarser levels. 00631 if (DoFixup) 00632 Fixup_().SetMaxRecomposeLevel(Level); 00633 register int lev; 00634 for (lev=Level; lev<=MaxLevel(GH()); lev++) { 00635 U().GF_DeleteStorage(1,lev); 00636 if (ShadowAllowed) { 00637 Ush().GF_DeleteStorage(1,lev); 00638 if (lev==MaxLevel(GH())) Ush().GF_DeleteStorage(0,lev); 00639 } 00640 } 00641 00642 #ifdef DEBUG_AMR 00643 if (Integrator_().Abort()) 00644 for (register int l=0; l<=FineLevel(GH()); l++) { 00645 int time = CurrentTime(GH(),l); 00646 int tstep = TimeStep(U(),l); 00647 CheckLevel(U(),time,l,"AMRSolver(DIM).AdvanceLevel::before Recompose"); 00648 if (l<Level) 00649 CheckLevel(U(),time+tstep,l, 00650 "AMRSolver(DIM).AdvanceLevel::before Recompose for Time+TStep"); 00651 if (ShadowAllowed && l<MaxLevel(GH())) { 00652 int shtstep = TimeStep(Ush(),l); 00653 CheckLevel(Ush(),time,l, 00654 "AMRSolver(DIM).AdvanceLevel::before Recompose for Shadow"); 00655 if (l<Level) 00656 CheckLevel(Ush(),time+shtstep,l, 00657 "AMRSolver(DIM).AdvanceLevel::before Recompose for Time+TStep for Shadow"); 00658 } 00659 } 00660 #endif 00661 00662 RecomposeHierarchy(GH(), (Level==0 && RecomposeBaseLev) || 00663 (Level>0 && RecomposeHighLev) ? DAGHTrue : DAGHFalse); 00664 #ifdef DEBUG_AMR 00665 if (Integrator_().Abort()) 00666 for (register int l=0; l<=FineLevel(GH()); l++) { 00667 int time = CurrentTime(GH(),l); 00668 int tstep = TimeStep(U(),l); 00669 CheckLevel(U(),time,l,"AMRSolver(DIM).AdvanceLevel::after Recompose"); 00670 if (l<Level) 00671 CheckLevel(U(),time+tstep,l, 00672 "AMRSolver(DIM).AdvanceLevel::after Recompose for Time+TStep"); 00673 if (ShadowAllowed && l<MaxLevel(GH())) { 00674 int shtstep = TimeStep(Ush(),l); 00675 CheckLevel(Ush(),time,l, 00676 "AMRSolver(DIM).AdvanceLevel::after Recompose for Shadow"); 00677 if (l<Level) 00678 CheckLevel(Ush(),time+shtstep,l, 00679 "AMRSolver(DIM).AdvanceLevel::after Recompose for Time+TStep for Shadow"); 00680 } 00681 } 00682 #endif 00683 00684 for (lev=Level; lev<=MaxLevel(GH()); lev++) { 00685 U().GF_CreateStorage(1,lev); 00686 if (ShadowAllowed && lev<MaxLevel(GH())) 00687 Ush().GF_CreateStorage(1,lev); 00688 } 00689 END_WATCH_RECOMPOSING_WHOLE 00690 00691 if (finest < FineLevel(GH())) { 00692 t[FineLevel(GH())] = t[finest]; 00693 dt[FineLevel(GH())] = dt[finest] / RefineFactor(GH(),finest); 00694 // Use this for simplified AMR 00695 // dt[FineLevel(GH())] = dt[finest]; 00696 } 00697 RegridDone = true; 00698 } 00699 else 00700 RegridDone = false; 00701 00702 //--------------------------------------------------------------- 00703 // Take one step on U(). U()(Time, Level) remains unchanged. 00704 // In principle, it would be possible to skip this step after error estimation on the 00705 // unchanged level. But this is not the case here, because numerical fluxes are only 00706 // available for single grids. As fine grids might change due to error estimation, 00707 // the saving of coarse grid fluxes would require the storage of numerical fluxes 00708 // on the entire coarse level. 00709 //--------------------------------------------------------------- 00710 cfl = IntegrateLevel(U(), Time, Level, t[Level], dt[Level], 00711 DoFixup, t[Level-1], 0); 00712 cfl_new[Level] = cfl_new[Level] > cfl ? cfl_new[Level] : cfl; 00713 00714 //--------------------------------------------------------------- 00715 // Take one step on Ush(). 00716 // If restriction into Ush() is done AFTER the step on 00717 // U(), U()(Time,Level) has to remain unchanged in a fractional 00718 // step method! 00719 //--------------------------------------------------------------- 00720 00721 if (RegridEvery != AMRSolverRegridFalse ? 00722 Level<MaxLevel(GH()) && ShadowAllowed && 00723 (Level>0 ? (k+1)%RegridEvery==0 : true) : false) { 00724 00725 double dt_sh = dt[Level] * Ush().factor(); 00726 SetPhysicalTime(Ush(),Time,Level,t[Level]); 00727 00728 //------------------------------------------------------------------------ 00729 // Restrict from U() into Ush() 00730 //------------------------------------------------------------------------ 00731 int myargc = sizeof(int); int myargs = NEquations(); 00732 forall (U(),Time,Level,cm) 00733 vec_grid_data_type& u_main = U()(Time, Level, cm); 00734 forall (Ush(),Time,Level,cs) 00735 vec_grid_data_type& u_shadow = Ush()(Time, Level, cs); 00736 BBox bb = u_shadow.bbox()*u_main.bbox(); 00737 if (!bb.empty()) 00738 f_restrict_amr(FA(DIM,u_main),FA(DIM,u_shadow), 00739 BOUNDING_BOX(bb),(char *)&myargs, &myargc); 00740 end_forall 00741 end_forall 00742 00743 START_WATCH 00744 AdaptiveBoundaryUpdate(Ush(),Time,Level); 00745 #ifdef DEBUG_AMR 00746 CheckLevel(Ush(),Time,Level, 00747 "AMRSolver(DIM).AdvanceLevel::after AdaptiveBoundaryUpdate for Shadow"); 00748 #endif 00749 END_WATCH_BOUNDARIES_INTERPOLATION 00750 START_WATCH 00751 Sync(Ush(),Time,Level); 00752 #ifdef DEBUG_AMR 00753 CheckLevel(Ush(),Time,Level,"AMRSolver(DIM).AdvanceLevel::after Sync for Shadow"); 00754 #endif 00755 END_WATCH_BOUNDARIES_SYNC 00756 START_WATCH 00757 ExternalBoundaryUpdate(Ush(),Time,Level); 00758 #ifdef DEBUG_AMR 00759 CheckLevel(Ush(),Time,Level, 00760 "AMRSolver(DIM).AdvanceLevel::after ExternalBoundaryUpdate for Shadow"); 00761 #endif 00762 END_WATCH_BOUNDARIES_WHOLE 00763 00764 IntegrateLevel(Ush(), Time, Level, t[Level], dt_sh, false, 0.0, 2); 00765 00766 } // Ends if statement for Ush() 00767 00768 //--------------------------------------------------------------- 00769 // If we are not at the finest level then go to next level 00770 //--------------------------------------------------------------- 00771 00772 if (Level < FineLevel(GH())) { 00773 00774 // Boundary conditions have to be applied on U()(Time+TStep,Level) 00775 // to ensure correct time-space interpolation for Level+1. 00776 START_WATCH 00777 AdaptiveBoundaryUpdate(U(),Time+TStep,Level); 00778 #ifdef DEBUG_AMR 00779 CheckLevel(U(),Time+TStep,Level, 00780 "AMRSolver(DIM).AdvanceLevel::after AdaptiveBoundaryUpdate for Time+TStep"); 00781 #endif 00782 END_WATCH_BOUNDARIES_INTERPOLATION 00783 00784 START_WATCH 00785 Sync(U(),Time+TStep,Level); 00786 #ifdef DEBUG_AMR 00787 CheckLevel(U(),Time+TStep,Level, 00788 "AMRSolver(DIM).AdvanceLevel::after Sync for Time+TStep"); 00789 #endif 00790 END_WATCH_BOUNDARIES_SYNC 00791 00792 START_WATCH 00793 ExternalBoundaryUpdate(U(),Time+TStep,Level); 00794 #ifdef DEBUG_AMR 00795 CheckLevel(U(),Time+TStep,Level, 00796 "AMRSolver(DIM).AdvanceLevel::after ExternalBoundaryUpdate for Time+TStep"); 00797 #endif 00798 END_WATCH_BOUNDARIES_WHOLE 00799 00800 // Recursive Integration of next level in hierarchy 00801 AdvanceLevel(Level+1, RegridEvery, RegridDone, 00802 ShadowAllowed, DoFixup, RecomposeBaseLev, 00803 RecomposeHighLev); 00804 00805 #ifdef DEBUG_AMR 00806 CheckLevel(U(),Time+TStep,Level,"AMRSolver(DIM).AdvanceLevel::before Restrict and Fixup"); 00807 #endif 00808 00809 #ifdef DEBUG_PRINT_FIXUP 00810 VectorType mass_old; 00811 if (Level == 0 && _Fixup) 00812 mass_old = Sum(U(),Time+TStep,Level); 00813 #endif 00814 00815 if (DoFixup) { 00816 START_WATCH 00817 Fixup_().Correction(Time+TStep, Time, Level, dx); 00818 END_WATCH_FIXUP_WHOLE 00819 } 00820 00821 // Implementation of fixup assumes that restriction is done AFTERWARDS! 00822 // Only values from internal cells are used. Ghostcell values are not restricted. 00823 int myargc = sizeof(int); int myargs = NEquations(); 00824 Restrict(U(), Time+TStep, Level+1, Time+TStep, Level, (char *)&myargs, myargc); 00825 00826 #ifdef DEBUG_AMR 00827 CheckLevel(U(),Time+TStep,Level,"AMRSolver(DIM).AdvanceLevel::after Restrict and Fixup"); 00828 #endif 00829 00830 #ifdef DEBUG_PRINT_FIXUP 00831 VectorType mass_new; 00832 if (Level == 0 && _Fixup) { 00833 ( comm_service::log() 00834 << "\n*********** Difference in state variables at level " 00835 << Level << " after fixup ********\n" ).flush(); 00836 // Sum() uses only internal cells. No synchronization on Level is therefore 00837 // necessary after restriction! 00838 mass_new = Sum(U(),Time+TStep,Level); 00839 ( comm_service::log() << mass_old << "-" 00840 << mass_new << "=" << mass_old - mass_new ).flush(); 00841 ( comm_service::log() 00842 << "\n**************************************************************" 00843 << "***********\n" ).flush(); 00844 } 00845 #endif 00846 00847 } 00848 #ifdef DEBUG_PROFILE 00849 ( cout << "\n" ).flush(); 00850 #endif 00851 00852 //--------------------------------------------------------------- 00853 // Move from current time to previous. 00854 //--------------------------------------------------------------- 00855 CycleTimeLevels(U(),Level); 00856 if (ShadowAllowed) 00857 CycleTimeLevels(Ush(),Level); 00858 00859 //--------------------------------------------------------------- 00860 // Increment time step on current level 00861 // IncrCurrentTime() moves internal time-indices. 00862 // Can only be called after all time-depended operations, like 00863 // time-interpolation or time-cycling. 00864 //--------------------------------------------------------------- 00865 t[Level] += dt[Level]; 00866 IncrCurrentTime(GH(),Level); 00867 SetPhysicalTime(U(),Time+TStep,Level,t[Level]); 00868 if (ShadowAllowed) 00869 SetPhysicalTime(Ush(),Time+TStep,Level,t[Level]+ 00870 dt[Level]*(Ush().factor()-1)); 00871 // Use this for simplified AMR 00872 // if (Level>0) 00873 // for (int kk=1; kk<RefineFactor(GH(),Level-1); kk++) 00874 // IncrCurrentTime(GH(),Level); 00875 RegridDone = false; 00876 } 00877 } 00878 00879 void CheckLevel(vec_grid_fct_type &u, const int time, const int level, const char *text) { 00880 if (Integrator_().Abort()) { 00881 char wtext[1024]; 00882 sprintf(wtext,"\non Level=%d at Time=%d: %s\n",level,time,text); 00883 forall (u,time,level,c) 00884 vec_grid_data_type& uw = u(time,level,c); 00885 Integrator_().CheckGrid(uw, uw.bbox(), wtext); 00886 end_forall 00887 } 00888 } 00889 00890 protected: 00891 BBox BoundaryBBox(const BBox& whole, const int s) { 00892 int d = s/2; 00893 #ifdef DEBUG_PRINT 00894 assert (d>=0 && d<DIM); 00895 #endif 00896 BBox bb(whole); 00897 if (s%2 == 0) 00898 bb.upper(d) = bb.lower(d)+(NGhosts()-1)*bb.stepsize(d); 00899 else 00900 bb.lower(d) = bb.upper(d)-(NGhosts()-1)*bb.stepsize(d); 00901 return bb; 00902 } 00903 00904 00905 int FixupPar; 00906 int RegridEvery; 00907 double MinEfficiency; 00908 int BlockWidth, OverlapWidth, BufferWidth, NestingBuffer; 00909 int PlotFlags; 00910 00911 double *t, *dt; 00912 double *cfl_new; 00913 }; 00914 00915 00916 #endif
Quickstart Users Guide Programmers Reference Installation Examples Download
AMROC Main Home Contactlast update: 06/01/04