amroc/amr/AMRSolverBase.h

Go to the documentation of this file.
00001 // -*- C++ -*-
00002 
00003 // Copyright (C) 2002 Ralf Deiterding
00004 // Brandenburgische Universitaet Cottbus
00005 //
00006 // Copyright (C) 2003-2007 California Institute of Technology
00007 // Ralf Deiterding, ralf@cacr.caltech.edu
00008 
00009 #ifndef AMROC_SOLVERBASE_H
00010 #define AMROC_SOLVERBASE_H
00011 
00019 #include "DAGH.h"
00020 #include "DAGHIO.h"             
00021 #include "DAGHCluster.h"
00022 
00023 #include "Solver.h"
00024 #include "Timing.h" 
00025 #include "Integrator.h"
00026 #include "InitialCondition.h"
00027 #include "BoundaryConditions.h"
00028 #include "LevelTransfer.h"
00029 #include "AMRFlagging.h"
00030 #include "AMRFixup.h"
00031 #include "ExactSolution.h"
00032 #include "FileOutput.h"
00033 #include "VectorGridDataOps.h"
00034 
00035 #include <iostream>
00036 #include <string>
00037 #include <cstdio>
00038 #include <cstdlib>
00039 #include <cfloat>
00040 #include <stdio.h>
00041 
00042 #define FixedTimeStepsNoCFL          (-1)
00043 #define FixedTimeSteps               (0)
00044 #define AutomaticTimeStepsNoRestart  (1)
00045 #define AutomaticTimeStepsRestart    (2)
00046 #define MaxCutOffRegions            (10)
00047 #define MaxAMRTimeSteps             (10)
00048 #define FACTOR 1.0e5
00049 
00050 class AMRTimeStep : public controlable {
00051 public:
00052   AMRTimeStep() : LastTime(1.e37), VariableTimeStepping(2) {
00053     dtv[0] = 0.01; 
00054     dtv[1] = 0.02;
00055     cflv[0] = 0.9;
00056     cflv[1] = 0.8;
00057   }
00058 
00059   virtual ~AMRTimeStep() {}
00060 
00061   virtual void register_at(ControlDevice& Ctrl) {}
00062   virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00063     LocCtrl = Ctrl.getSubDevice(prefix+"TimeControl");
00064     RegisterAt(LocCtrl,"LastTime",LastTime);
00065     RegisterAt(LocCtrl,"StepMode",VariableTimeStepping);
00066     RegisterAt(LocCtrl,"TimeStep",dtv[0]);
00067     RegisterAt(LocCtrl,"TimeStepMax",dtv[1]);
00068     RegisterAt(LocCtrl,"CFLRestart",cflv[0]);
00069     RegisterAt(LocCtrl,"CFLControl",cflv[1]);
00070   }
00071   
00072 public: 
00073   double LastTime;
00074   int VariableTimeStepping;
00075   double dtv[2], cflv[2];
00076   ControlDevice LocCtrl;     
00077 };
00078 
00079 
00089 template <class VectorType, class FixupType, class FlagType, int dim>
00090 class AMRSolverBase : public Solver {
00091   typedef AMRBase<VectorType,dim> base;
00092   typedef typename VectorType::InternalDataType DataType;
00093 public:
00094   typedef GridFunction<VectorType,dim> vec_grid_fct_type;  
00095   typedef GridData<VectorType,dim> vec_grid_data_type;  
00096   typedef GridFunction<DataType,dim> grid_fct_type;
00097   typedef GridData<DataType,dim> grid_data_type;
00098 
00099   typedef Integrator<VectorType,dim> integrator_type;
00100   typedef InitialCondition<VectorType,dim> initial_condition_type;
00101   typedef BoundaryConditions<VectorType,dim> boundary_conditions_type;
00102   typedef LevelTransfer<VectorType,dim> leveltransfer_type;
00103   typedef AMRFlagging<VectorType,FixupType,FlagType,dim> flagging_type;
00104   typedef AMRFixup<VectorType,FixupType,dim> fixup_type;
00105   typedef ExactSolution<VectorType,dim> exact_solution_type;
00106   typedef FileOutput<VectorType,dim> file_output_type;
00107 
00108   typedef GFBndryUpdateSpecificFunc<boundary_conditions_type,VectorType,dim> boundary_functor_type;
00109   typedef GFLevelTransferSpecificFunc<leveltransfer_type,VectorType,dim> leveltransfer_functor_type;
00110   typedef GFAdptBndryUpdateSpecificFunc<leveltransfer_type,VectorType,dim> adaptbnd_functor_type;
00111 
00112   AMRSolverBase(integrator_type& integ, initial_condition_type& init, boundary_conditions_type& bc) : 
00113     _Integrator(integ), _InitialCondition(init), _BoundaryConditions(bc),
00114     _Equations(VectorType::Length()), _Ghosts(2), _Dim(dim), RedistributeEvery(0), 
00115     Distribution(DAGHCompositeDistribution), MaxLev(1), GuCFactor(2), MaxGridBoxSize(0), CutOffs(0), 
00116     NAMRTimeSteps(1), firstnode(0), lastnode(-1), UseIOServer(0) {
00117     for (register int d=0; d<dim; d++) {
00118       shape[d] = 2;
00119       geom[2*d] = 0.0;
00120       geom[2*d+1] = 1.0;
00121       periodic[d] = DAGHFalse;
00122     }    
00123     for (register int i=0; i<DAGHMaxLevels; i++)
00124       _RefineFactor[i] = DAGHDefaultRefineFactor;
00125 
00126     _Hierarchy = (GridHierarchy*) 0;
00127     _u = (vec_grid_fct_type *) 0;
00128     _u_sh = (vec_grid_fct_type *) 0;
00129     _work = (grid_fct_type *) 0;
00130     _work_sh = (grid_fct_type *) 0;
00131     _LevelTransfer = (leveltransfer_type *) 0;
00132     _Flagging = (flagging_type *) 0;
00133     _Fixup = (fixup_type *) 0;
00134     _ExactSolution = (exact_solution_type *) 0;
00135     _FileOutput = (file_output_type *) 0;
00136     _BoundaryFunc = (boundary_functor_type *) 0;
00137     _ProlongFunc = (leveltransfer_functor_type *) 0;
00138     _RestrictFunc = (leveltransfer_functor_type *) 0;
00139     _AdaptBndFunc = (adaptbnd_functor_type *) 0;
00140     _LastOutputTime = -1;
00141     CheckpointName = "checkpt";
00142     CheckpointSave = 1;
00143     FinalWaitingTime = 10.;
00144   }
00145 
00146   virtual ~AMRSolverBase() {
00147     if (_u) delete _u;
00148     if (_u_sh) delete _u_sh;
00149     if (_work) delete _work;
00150     if (_work_sh) delete _work_sh;
00151     if (_BoundaryFunc) delete _BoundaryFunc;
00152     if (_ProlongFunc) delete _ProlongFunc;
00153     if (_RestrictFunc) delete _RestrictFunc;
00154     if (_AdaptBndFunc) delete _AdaptBndFunc;
00155     if (_Hierarchy) delete _Hierarchy;
00156   }
00157 
00158   //******************************************************************************
00159   // Abstract class interface
00160   //******************************************************************************
00161   virtual double Tick(int VariableTimeStepping, const double dtv[], const double cflv[],
00162                       int& Rejections) = 0;
00163   virtual double IntegrateLevel(vec_grid_fct_type& u, const int Time, const int Level, 
00164                                 double t, double dt, bool DoFixup, double tc, 
00165                                 const int which) = 0;
00166   virtual void Initialize_(const double& dt_start) = 0;    
00167   virtual bool Restart_(const char* CheckpointFile) = 0;  
00168   virtual void Checkpointing_(const char* CheckpointFile) = 0;  
00169   //******************************************************************************
00170 
00171   virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00172     PCtrl = Ctrl.getSubDevice("LocalGroup");
00173     RegisterAt(PCtrl,"FirstNode",firstnode);
00174     RegisterAt(PCtrl,"LastNode",lastnode);
00175     RegisterAt(PCtrl,"UseIOServer",UseIOServer);
00176 
00177     GHCtrl = Ctrl.getSubDevice("GridHierarchy");
00178     char VariableName[32];
00179     for (int d=0; d<dim; d++) {
00180       std::sprintf(VariableName,"Cells(%d)",d+1);
00181       RegisterAt(GHCtrl,VariableName,shape[d]);
00182       std::sprintf(VariableName,"GeomMin(%d)",d+1);
00183       RegisterAt(GHCtrl,VariableName,geom[2*d]);
00184       std::sprintf(VariableName,"GeomMax(%d)",d+1);
00185       RegisterAt(GHCtrl,VariableName,geom[2*d+1]);     
00186       std::sprintf(VariableName,"PeriodicBoundary(%d)",d+1);
00187       RegisterAt(GHCtrl,VariableName,periodic[d]);     
00188     }
00189     for (int nc=0; nc<MaxCutOffRegions; nc++) {
00190       for (int d=0; d<dim; d++) 
00191         for (int j=0; j<2; j++) {
00192           std::sprintf(VariableName,"CCells(%d,%d,%d)",nc+1,j+1,d+1);
00193           RegisterAt(GHCtrl,VariableName,cuts[nc*2*dim+2*d+j]); 
00194           cuts[nc*2*dim+2*d+j] = -100000;
00195         }
00196     }
00197     RegisterAt(GHCtrl,"CutOffs",CutOffs); 
00198     for (register int i=0; i<DAGHMaxLevels; i++) {
00199       std::sprintf(VariableName,"RefineFactor(%d)",i+1);
00200       RegisterAt(GHCtrl,VariableName,_RefineFactor[i]);
00201     }
00202     RegisterAt(GHCtrl,"RedistributeEvery",RedistributeEvery);
00203     RegisterAt(GHCtrl,"MaxLevels",MaxLev);
00204     RegisterAt(GHCtrl,"Distribution",Distribution);
00205     RegisterAt(GHCtrl,"GuCFactor",GuCFactor); 
00206     RegisterAt(GHCtrl,"MaxGridBoxSize",MaxGridBoxSize); 
00207     RegisterAt(GHCtrl,"GhostCells",_Ghosts);
00208     RegisterAt(GHCtrl,"CheckpointName",CheckpointName);    
00209     RegisterAt(GHCtrl,"CheckpointSave",CheckpointSave);
00210     RegisterAt(GHCtrl,"FinalWaitingTime",FinalWaitingTime);
00211   
00212     Step[0].register_at(Ctrl,"");
00213     for (int cs=0; cs<MaxAMRTimeSteps; cs++) {
00214       char text[5];
00215       std::sprintf(text,"%d-",cs+1);
00216       Step[cs].register_at(Ctrl,std::string(text));
00217     }
00218     RegisterAt(Ctrl,"TimeControls",NAMRTimeSteps); 
00219 
00220     _Integrator.register_at(Ctrl, "");
00221     _InitialCondition.register_at(Ctrl, "");
00222     _BoundaryConditions.register_at(Ctrl, ""); 
00223     if (_LevelTransfer) _LevelTransfer->register_at(Ctrl, ""); 
00224     if (_Flagging) _Flagging->register_at(Ctrl, ""); 
00225     if (_Fixup) _Fixup->register_at(Ctrl, ""); 
00226     if (_ExactSolution) _ExactSolution->register_at(Ctrl, ""); 
00227     if (_FileOutput) _FileOutput->register_at(Ctrl, ""); 
00228   }
00229   virtual void register_at(ControlDevice& Ctrl) { register_at(Ctrl, ""); }
00230 
00231   virtual void init() {
00232     _Integrator.init();
00233     _InitialCondition.init();
00234     _BoundaryConditions.init();
00235     if (_LevelTransfer) _LevelTransfer->init();
00236     if (_Flagging) _Flagging->init();
00237     if (_Fixup) _Fixup->init();
00238     if (_ExactSolution) _ExactSolution->init();
00239     if (_FileOutput) _FileOutput->init();
00240   }
00241   virtual void update() {
00242     _Integrator.update();
00243     _InitialCondition.update();
00244     _BoundaryConditions.update();
00245     if (_LevelTransfer) _LevelTransfer->update();
00246     if (_Flagging) _Flagging->update();
00247     if (_Fixup) _Fixup->update();
00248     if (_ExactSolution) _ExactSolution->update();
00249     if (_FileOutput) _FileOutput->update();
00250     START_WATCH_WHOLE
00251   }
00252   virtual void finish() {
00253     END_WATCH_WHOLE
00254 #ifdef TIMING_AMR
00255     if (_Hierarchy) {
00256       Timing::collect(Comm());
00257 #ifdef DEBUG_PRINT
00258       Timing::print_local(comm_service::log());
00259 #endif
00260       if (MY_PROC == VizServer)
00261         Timing::print(std::cout);
00262     }
00263 #endif
00264     _Integrator.finish();
00265     _InitialCondition.finish();
00266     _BoundaryConditions.finish();
00267     if (_LevelTransfer) _LevelTransfer->finish();
00268     if (_Flagging) _Flagging->finish();
00269     if (_Fixup) _Fixup->finish();
00270     if (_ExactSolution) _ExactSolution->finish();
00271     if (_FileOutput) _FileOutput->finish();
00272     if (_u) {
00273       delete _u;
00274       _u = (vec_grid_fct_type *) 0;
00275     }
00276     if (_u_sh) { 
00277       delete _u_sh;
00278       _u_sh = (vec_grid_fct_type *) 0;
00279     }
00280     if (_work) {
00281       delete _work;
00282       _work = (grid_fct_type *) 0;
00283     }
00284     if (_work_sh) {
00285       delete _work_sh;
00286       _work_sh = (grid_fct_type *) 0;
00287     }
00288     if (_BoundaryFunc) {
00289       delete _BoundaryFunc;
00290       _BoundaryFunc = (boundary_functor_type *) 0;
00291     }
00292     if (_ProlongFunc) {
00293       delete _ProlongFunc;
00294       _ProlongFunc = (leveltransfer_functor_type *) 0;
00295     }
00296     if (_RestrictFunc) {
00297       delete _RestrictFunc;
00298       _RestrictFunc = (leveltransfer_functor_type *) 0;
00299     }
00300     if (_AdaptBndFunc) {
00301       delete _AdaptBndFunc;
00302       _AdaptBndFunc = (adaptbnd_functor_type *) 0;
00303     }
00304 #ifndef DAGH_NO_MPI
00305     if (FinalWaitingTime>0. && (_Hierarchy!=(GridHierarchy *) 0)) {
00306       double start_finalize = MPI_Wtime();
00307       while (MPI_Wtime()-start_finalize < FinalWaitingTime) {}
00308     }
00309 #endif
00310 
00311     if (_Hierarchy) delete _Hierarchy;
00312     _Hierarchy = (GridHierarchy *) 0;
00313 
00314 #ifndef DAGH_NO_MPI
00315     comm_service::clean();
00316 #endif
00317   } 
00318 
00319   virtual bool setup() {
00320     bool ret = true;
00321     //---------------------------------------------------------------
00322     // Reserve last node for IO
00323     //---------------------------------------------------------------
00324     if (UseIOServer) comm_service::set_io_enable();
00325     else comm_service::reset_io_enable();
00326 
00327 #ifndef DAGH_NO_MPI
00328     int size, rank = 0;
00329     MPI_Comm_size(MPI_COMM_WORLD, &size);
00330     MPI_Comm_rank(MPI_COMM_WORLD, &rank);
00331 
00332     if (firstnode<0) firstnode=0;
00333     if (lastnode>size-1 || lastnode<0) lastnode=size-1;
00334 
00335     if (firstnode<=rank && rank<=lastnode) {
00336       MPI_Group GrpWorld, GrpLocal;
00337       MPI_Comm CommLocal;
00338       MPI_Comm_group(MPI_COMM_WORLD, &GrpWorld);
00339       
00340       int* localranks = new int [lastnode-firstnode+1];
00341       int localcnt = 0;
00342       for (register int i=firstnode; i<=lastnode; i++)
00343         localranks[localcnt++] = i;
00344       
00345       MPI_Group_incl(GrpWorld, localcnt, localranks, &GrpLocal);
00346       MPI_Comm_create(MPI_COMM_WORLD, GrpLocal, &CommLocal);
00347       
00348       comm_service::init(CommLocal);
00349 
00350       SetupData();
00351       ret = COMPUTE_NODE;
00352     }
00353     else 
00354       ret = false;
00355 #else
00356     firstnode=0;
00357     lastnode=0;    
00358     SetupData();
00359     ret = COMPUTE_NODE;
00360 #endif
00361     return ret;
00362   }
00363 
00364   virtual void SetupGridHierarchy() {
00365     _Hierarchy = new GridHierarchy(dim,DAGHCellCentered,MaxLev);
00366     SetBaseGrid(GH(),geom,shape,CutOffs,cuts);   
00367     for (register int l=0; l<MaxLev-1; l++) {
00368       if (_RefineFactor[l]<DAGHDefaultRefineFactor)
00369         _RefineFactor[l] = DAGHDefaultRefineFactor;
00370       if (ErrorEstimation() && _RefineFactor[l]%DAGHShadowFactor && l<MaxLev-2) {
00371         std::cout << "   RefineFactor(" << l+1
00372                   << ") is not dividable by Shadowfactor()." 
00373                   << " Can't use error estimation.\n" << std::flush;
00374         Flagging_().SetErrorEstimation(false);
00375       } 
00376       SetRefineFactor(GH(),l,_RefineFactor[l]); 
00377     }
00378 
00379     SetBoundaryWidth(GH(),NGhosts());
00380     SetBoundaryType(GH(),DAGHBoundaryUserDef); 
00381     if (ErrorEstimation()) 
00382       GuCFactor = Max(GuCFactor,DAGHShadowFactor);
00383     if (GuCFactor) SetGuCFactor(GH(),GuCFactor);
00384     if (MaxGridBoxSize) SetMaxGridBoxSize(GH(), MaxGridBoxSize);
00385     for (int d=0; d<dim; d++)
00386       if (periodic[d]) SetPeriodicBoundaries(GH(),d,DAGHTrue);
00387     SetDistributionType(GH(),Distribution);
00388     
00389     if (_FileOutput) _FileOutput->SetupData(PGH(), NGhosts());
00390 
00391     ComposeHierarchy(GH());     
00392     _Integrator.SetupData(PGH(), NGhosts());
00393     _InitialCondition.SetupData(PGH(), NGhosts());
00394     _BoundaryConditions.SetupData(PGH(), NGhosts());
00395     if (_LevelTransfer) _LevelTransfer->SetupData(PGH(), NGhosts());
00396     if (_Flagging) _Flagging->SetupData(PGH(), NGhosts());
00397     if (_Fixup) _Fixup->SetupData(PGH(), NGhosts());
00398     if (_ExactSolution) _ExactSolution->SetupData(PGH(), NGhosts());
00399   }
00400 
00401   virtual void SetupData() {
00402     START_WATCH;
00403     if (_Hierarchy==(GridHierarchy *) 0) 
00404       SetupGridHierarchy();
00405 
00406     //---------------------------------------------------------------
00407     //  Declare and register available functors for GridFunctions
00408     //---------------------------------------------------------------
00409     _BoundaryFunc = new boundary_functor_type(&_BoundaryConditions, &boundary_conditions_type::SetBndry);
00410     if (_LevelTransfer) {
00411       _ProlongFunc = new leveltransfer_functor_type(_LevelTransfer, &leveltransfer_type::Prolong);
00412       _RestrictFunc = new leveltransfer_functor_type(_LevelTransfer, &leveltransfer_type::Restrict);
00413       if (_LevelTransfer->UseAdaptBndry()) 
00414         _AdaptBndFunc = new adaptbnd_functor_type(_LevelTransfer, &leveltransfer_type::SetAdaptBndry);
00415     }
00416 
00417     //---------------------------------------------------------------
00418     //  Declare GridFunctions    
00419     //---------------------------------------------------------------
00420     AllocError::SetTexts("Solver::SetupData","allocation of GridFunctions");
00421     int t_sten = 1;   
00422     int s_sten = NGhosts();       
00423     std::sprintf(GFName,"u");
00424     _u = new vec_grid_fct_type(GFName, t_sten, s_sten, GH(), DAGHCommSimple);
00425     SetTimeAlias(U(), -1, 1);
00426     SetBndryUpdateFunction(U(), _BoundaryFunc);
00427     if (_LevelTransfer) {
00428       SetProlongFunction(U(), _ProlongFunc);
00429       SetRestrictFunction(U(), _RestrictFunc);
00430       if (_LevelTransfer->UseAdaptBndry()) {
00431         U().GF_SetAdaptBoundaryType(_LevelTransfer->AdaptiveBoundaryType());
00432         U().GF_SetAdaptiveBndryUpdateFunc(_AdaptBndFunc);
00433       }
00434     }
00435 
00436     if (ErrorEstimation()) {
00437       std::sprintf(GFNamesh,"ush");     
00438       _u_sh = new vec_grid_fct_type(GFNamesh, t_sten, s_sten,  
00439                                     DAGHShadowFactor, GH(), DAGHCommSimple);
00440       SetTimeAlias(Ush(), -1, 1);
00441       SetBndryUpdateFunction(Ush(), _BoundaryFunc);
00442       if (_LevelTransfer) {
00443         SetProlongFunction(Ush(), _ProlongFunc);
00444         SetRestrictFunction(Ush(), _RestrictFunc);
00445         if (_LevelTransfer->UseAdaptBndry()) {
00446           Ush().GF_SetAdaptBoundaryType(_LevelTransfer->AdaptiveBoundaryType());
00447           Ush().GF_SetAdaptiveBndryUpdateFunc(_AdaptBndFunc);
00448         }
00449       }
00450     }
00451       
00452     // One work array with two time steps and Shadow-hierarchy is typically necessary 
00453     std::sprintf(IOName,"IOPlaceholder"); 
00454     _work = new grid_fct_type(IOName, t_sten, s_sten, GH(), DAGHCellCentered, DAGHCommSimple, 
00455                              DAGHNoBoundary, DAGHNoAdaptBoundary, DAGHNoExternalGhost);
00456     SetCheckpointFlag(Work(), DAGHFalse);
00457     Work().GF_SetMaxRecomposeLevel(DAGHNull);
00458     SetTimeAlias(Work(), -1, 1);
00459  
00460     if (ErrorEstimation()) {
00461       std::sprintf(IONamesh,"IOPlaceholdersh"); 
00462       _work_sh = new grid_fct_type(IONamesh, t_sten, s_sten, GH(), DAGHCellCentered, 0, 
00463                                    DAGHShadowFactor, DAGH_All, DAGHCommSimple, 
00464                                    DAGHNoBoundary, DAGHNoAdaptBoundary, DAGHNoExternalGhost);
00465       SetCheckpointFlag(Worksh(), DAGHFalse);
00466       Worksh().GF_SetMaxRecomposeLevel(DAGHNull);
00467       SetTimeAlias(Worksh(), -1, 1);
00468     }
00469     END_WATCH_COMPOSING_WHOLE
00470   } 
00471 
00472   virtual void Initialize(double& t, double &dt) {
00473     int me = MY_PROC;
00474     if (me == VizServer)    
00475       std::cout << " --- Initializing --- " << std::flush;
00476     Initialize_(Step[0].dtv[0]);  
00477     if (me == VizServer)    
00478       std::cout << " Levels=" << FineLevel(GH())+1 << std::endl;
00479     double Current_Time, Next_Time;
00480     GH().DAGH_GetTimeSpecs(Current_Time, Next_Time);
00481     if (_ExactSolution) 
00482       _ExactSolution->ErrorNorm(U(), Work(), Current_Time);
00483     t = Current_Time;
00484     dt = Next_Time-Current_Time;
00485   }
00486 
00487   virtual void Restart(double& t, double &dt) {
00488     int me = MY_PROC;
00489     if (me == VizServer)     
00490       std::cout << " --- Restarting from " << CheckpointName.c_str()
00491                 << " --- " << std::flush;
00492     if (Restart_(CheckpointName.c_str())) {
00493       if (me == VizServer)    
00494         std::cout << "   Last Iteration: " << StepsTaken(GH(),0) << std::endl;
00495 
00496       double Current_Time, Next_Time;
00497       GH().DAGH_GetTimeSpecs(Current_Time, Next_Time);
00498       if (_ExactSolution) 
00499         _ExactSolution->ErrorNorm(U(), Work(), Current_Time);
00500       t = Current_Time;
00501       dt = Next_Time-Current_Time;
00502     }
00503     else {
00504       if (me == VizServer)    
00505         std::cout << "   File not found. " << std::endl;
00506     }
00507   }
00508 
00509   virtual void Advance(double& t, double& dt) {
00510     int me = MY_PROC;
00511     double Current_Time, Next_Time, Final_Time;
00512     GH().DAGH_GetTimeSpecs(Current_Time, Next_Time);
00513     Final_Time = Current_Time + dt;
00514 
00515     int CStepControl = 0;    
00516 
00517     while ((Final_Time-Step[CStepControl].LastTime)/Final_Time > DBL_EPSILON*FACTOR && CStepControl<NAMRTimeSteps)
00518       CStepControl++;  
00519     if (CStepControl>=NAMRTimeSteps) CStepControl = NAMRTimeSteps-1;
00520 
00521     if ((Step[CStepControl].LastTime-Current_Time)/Step[CStepControl].LastTime > DBL_EPSILON*FACTOR ||
00522         CStepControl==NAMRTimeSteps-1) {   
00523 
00524       if (Next_Time > Final_Time)
00525         SetTimeSpecs(GH(), Current_Time, Final_Time, 2);
00526       
00527       if (me == VizServer) { 
00528         std::cout << " --- Iteration: " << StepsTaken(GH(),0)+1 << " ---";
00529         std::cout << "   t = " << Current_Time << std::flush;  
00530         } 
00531       
00532       int Rejections;
00533       double cfl_max = Tick(Step[CStepControl].VariableTimeStepping, 
00534                             Step[CStepControl].dtv, Step[CStepControl].cflv,
00535                             Rejections);
00536 
00537       double Help_Time;
00538       GH().DAGH_GetTimeSpecs(Next_Time, Help_Time);
00539       if (Help_Time > Step[CStepControl].LastTime) {   
00540         Help_Time = Step[CStepControl].LastTime;
00541       }
00542       if (me == VizServer) { 
00543         std::cout << "  dt=" << Next_Time-Current_Time 
00544                   << "  Levels=" << FineLevel(GH())+1;
00545         if (Step[CStepControl].VariableTimeStepping >= 0)
00546           std::cout << "  max. CFL=" << cfl_max;
00547       }
00548       Current_Time = Next_Time;
00549       Next_Time = Help_Time;
00550       
00551       if (me == VizServer) {
00552         if (Rejections)
00553           std::cout << "  Restarted: " << Rejections << "x";
00554         std::cout << std::endl << std::flush;
00555         }
00556       if (_ExactSolution) 
00557         _ExactSolution->ErrorNorm(U(), Work(), Current_Time);
00558     }  
00559 
00560     t = Current_Time;
00561     dt = Next_Time-Current_Time;
00562   }  
00563 
00564   virtual void Output() {
00565     if (CurrentTime(GH(),0) == LastOutputTime()) return;
00566     if (_FileOutput) {
00567       START_WATCH 
00568         _FileOutput->WriteOut(U(), Work());     
00569       END_WATCH_OUPUT  
00570       SetLastOutputTime(CurrentTime(GH(),0));
00571     }
00572     else 
00573       SetLastOutputTime(-1);
00574   }
00575 
00576   virtual void Checkpointing() {
00577     if (CurrentTime(GH(),0) == LastCheckpointTime()) return;
00578     int me = MY_PROC;
00579     if (me == VizServer) 
00580       std::cout << " *** Checkpoint at " << CurrentTime(GH(),0) 
00581                 << ": " << std::flush; 
00582     START_WATCH 
00583 #ifndef DAGH_NO_MPI
00584       MPI_Barrier(comm_service::comm()); 
00585 #endif
00586       if (_FileOutput) _FileOutput->CloseIO();
00587       if (CheckpointSave) {
00588         if (me == VizServer) 
00589           std::cout << "Moving..." << std::flush; 
00590         std::ostringstream newName_str;
00591         GH().chkpt_filename(newName_str,CheckpointName.c_str(),me);
00592         char oldName[256], newName[256];
00593         int len = newName_str.str().size();
00594         newName_str.str().copy(newName,len);
00595         int last = -1;
00596         for (int i=0; i<len-1; i++) 
00597           if (newName[i]=='/') 
00598             last = i;
00599         if (last>=0) std::strncpy(oldName,newName,last+1);
00600         oldName[++last] = '.';
00601         std::strncpy(oldName+last+1,newName+last,len-last);
00602         oldName[len+1] = '\0';
00603         rename(newName,oldName);
00604 #ifndef DAGH_NO_MPI
00605         MPI_Barrier(comm_service::comm()); 
00606 #endif
00607         if (me == VizServer) 
00608           std::cout << "done. " << std::flush; 
00609       }
00610       if (me == VizServer) 
00611         std::cout << "Writing..." << std::flush; 
00612       Checkpointing_(CheckpointName.c_str());
00613 #ifndef DAGH_NO_MPI
00614       MPI_Barrier(comm_service::comm()); 
00615 #endif
00616       if (me == VizServer) 
00617         std::cout << "done." << std::endl << std::flush; 
00618     END_WATCH_OUPUT  
00619     SetLastCheckpointTime(CurrentTime(GH(),0));
00620   }
00621 
00622   virtual int NSteps() { return StepsTaken(GH(),0); }
00623 
00624   GridHierarchy& GH() { return *_Hierarchy; }
00625   const GridHierarchy& GH() const { return *_Hierarchy; }
00626   inline GridHierarchy* PGH() const { return _Hierarchy; }
00627 
00628   const int& NEquations() const { return _Equations; }
00629   const int& NGhosts() const { return _Ghosts; }
00630   const int& Dim() const { return _Dim; }
00631 
00632   inline vec_grid_fct_type& U() { return *_u; }
00633   inline const vec_grid_fct_type& U() const { return *_u; }
00634   inline vec_grid_fct_type& Ush() { return *_u_sh; }
00635   inline const vec_grid_fct_type& Ush() const { return *_u_sh; }
00636   inline grid_fct_type& Work() { return *_work; }
00637   inline const grid_fct_type& Work() const { return *_work; }
00638   inline grid_fct_type& Worksh() { return *_work_sh; }
00639   inline const grid_fct_type& Worksh() const { return *_work_sh; }
00640 
00641   inline integrator_type& Integrator_() { return _Integrator; }
00642   inline const integrator_type& Integrator_() const { return _Integrator; }
00643 
00644   inline initial_condition_type& InitialCondition_() { return _InitialCondition; }
00645   inline const initial_condition_type& InitialCondition_() const 
00646     { return _InitialCondition; }
00647 
00648   inline boundary_conditions_type& BoundaryConditions_() { return _BoundaryConditions; }
00649   inline const boundary_conditions_type& BoundaryConditions_() const 
00650     { return _BoundaryConditions; }
00651 
00652   inline void SetLevelTransfer(leveltransfer_type* _leveltransfer) 
00653   { _LevelTransfer = _leveltransfer; }
00654   inline leveltransfer_type& LevelTransfer_() { return *_LevelTransfer; }
00655   inline const leveltransfer_type& LevelTransfer_() const { return *_LevelTransfer; }
00656 
00657   inline void SetFlagging(flagging_type* _flagging) { _Flagging = _flagging; }
00658   inline flagging_type& Flagging_() { return *_Flagging; }
00659   inline const flagging_type& Flagging_() const { return *_Flagging; }
00660   inline bool ErrorEstimation() const { 
00661     if (_Flagging)
00662       return Flagging_().ErrorEstimation();
00663     else return false;
00664   }
00665 
00666   inline void SetFixup(fixup_type* _fixup) { _Fixup = _fixup; }
00667   inline fixup_type& Fixup_() { return *_Fixup; }
00668   inline const fixup_type& Fixup_() const { return *_Fixup; }
00669   inline bool FixupPossible() const { return (_Fixup!=(fixup_type *)0); }
00670 
00671   inline void SetExactSolution(exact_solution_type* exact) { _ExactSolution = exact; }
00672   inline exact_solution_type& ExactSolution_() { return *_ExactSolution; }
00673   inline const exact_solution_type& ExactSolution_() const { return *_ExactSolution; }
00674 
00675   inline void SetFileOutput(file_output_type* output) { _FileOutput = output; }
00676   inline file_output_type& FileOutput_() { return *_FileOutput; }
00677   inline const file_output_type& FileOutput_() const { return *_FileOutput; }
00678 
00679   inline void SetLastOutputTime(int Time) { _LastOutputTime = Time; }
00680   inline int LastOutputTime() const { return _LastOutputTime; }
00681   inline void SetLastCheckpointTime(int Time) { _LastCheckpointTime = Time; }
00682   inline int LastCheckpointTime() const { return _LastCheckpointTime; }
00683 
00684   virtual int NMethodOrder() const { return Integrator_().NMethodOrder(); } 
00685 
00686   inline int FirstNode() const { return firstnode; }
00687   inline int LastNode() const { return lastnode; }
00688   inline int NNodes() const { return lastnode-firstnode+1; }
00689 
00690   inline int SizeWorld() const { return comm_service::proc_world(); }
00691   inline int Size() const { return comm_service::proc_num(); }
00692 
00693   inline MPI_Group GrpWorld() const { return comm_service::grp_world(); }
00694   inline MPI_Group Grp() const { return comm_service::grp(); }
00695 
00696   inline MPI_Comm CommWorld() const { return comm_service::comm_world(); }
00697   inline MPI_Comm Comm() const { return comm_service::comm(); }
00698 
00699 protected:
00700   integrator_type& _Integrator;
00701   initial_condition_type& _InitialCondition;
00702   boundary_conditions_type& _BoundaryConditions;
00703   int _Equations, _Ghosts, _Dim;
00704   GridHierarchy* _Hierarchy;
00705   vec_grid_fct_type *_u, *_u_sh;
00706   grid_fct_type *_work, *_work_sh;
00707   leveltransfer_type* _LevelTransfer;
00708   flagging_type* _Flagging;
00709   fixup_type* _Fixup;
00710   exact_solution_type* _ExactSolution;
00711   file_output_type* _FileOutput;
00712   boundary_functor_type* _BoundaryFunc;
00713   leveltransfer_functor_type *_ProlongFunc, *_RestrictFunc;  
00714   adaptbnd_functor_type* _AdaptBndFunc;
00715   int RedistributeEvery;
00716   int Distribution;  
00717   int MaxLev;   
00718   int GuCFactor;
00719   int MaxGridBoxSize;
00720   int _RefineFactor[DAGHMaxLevels];
00721   int periodic[dim];
00722   int shape[dim];  
00723   double geom[2*dim];
00724   int CutOffs;
00725   int cuts[2*dim*MaxCutOffRegions];
00726   int NAMRTimeSteps;            
00727   AMRTimeStep Step[MaxAMRTimeSteps];
00728   int _LastOutputTime, _LastCheckpointTime;
00729   int firstnode, lastnode, UseIOServer;
00730   std::string CheckpointName;
00731   int CheckpointSave;  
00732   double FinalWaitingTime;
00733   ControlDevice LocCtrl, GHCtrl, PCtrl;
00734   char GFName[DAGHBktGFNameWidth], GFNamesh[DAGHBktGFNameWidth];
00735   char IOName[DAGHBktGFNameWidth], IONamesh[DAGHBktGFNameWidth];
00736 };
00737 
00738 
00739 #endif
Generated on Thu Jun 30 02:19:13 2016 for AMROC Fluid-solver Framework by 1.6.3