AMROC - Source code documentation

00001 #ifndef AMROC_SOLVER_CONTROL_H
00002 #define AMROC_SOLVER_CONTROL_H
00003 
00012 #define MaxCutOffRegions 10
00013 #define MaxStepControls  10
00014 
00015 //---------------------------------------------------------------------- 
00016 // Flags              
00017 //---------------------------------------------------------------------- 
00018 #define VizServer 0                                               
00019 
00020 #include "iostream.h" 
00021 #include <float.h>
00022 
00023 #define FACTOR 1.0e5
00024 
00025 //---------------------------------------------------------------------- 
00026 // Timing includes             
00027 //---------------------------------------------------------------------- 
00028 #include "Timing.h" 
00029 #include "Timing.C"  
00030 
00031 //---------------------------------------------------------------------- 
00032 // DAGH includes             
00033 //---------------------------------------------------------------------- 
00034 #include "DAGH.h"
00035 #include "DAGHIO.h"             
00036 #include "DAGHCluster.h"
00037 
00038 #include "Solver.h"
00039 #include "ExactSolution.h"
00040 #include "SolverControlInterface.h"
00041 
00042 #ifdef SPX
00043 #   define _H_UNISTD
00044 #endif
00045 
00046 extern "C" {
00047 #include "sds.h"
00048 #include "amrsds.h"
00049 }
00050     
00051 #ifndef DAGH_NO_MPI 
00052 #include "mpi.h"
00053 #endif 
00054 
00055 #ifndef SolverControlName
00056 #define SolverControl(dim)      name2(SolverControl,dim)
00057 #define SolverControlName
00058 #endif
00059 
00060 class StepControl : public controlable {
00061 public:
00062   StepControl() : LastTime(0.0), VariableTimeStepping(0), Outputs(0) {
00063     dtv[0] = 0.01; 
00064     dtv[1] = 0.0;
00065     cflv[0] = 0.9;
00066     cflv[1] = 0.8;
00067   }
00068 
00069   virtual void register_at(ControlDevice& Ctrl) {}
00070   virtual void register_at(ControlDevice& Ctrl, const string& prefix) {
00071     LocCtrl = Ctrl.getSubDevice(prefix+"Time");
00072     RegisterAt(LocCtrl,"LastTime",LastTime);
00073     RegisterAt(LocCtrl,"StepMode",VariableTimeStepping);
00074     RegisterAt(LocCtrl,"TimeStep",dtv[0]);
00075     RegisterAt(LocCtrl,"TimeStepMax",dtv[1]);
00076     RegisterAt(LocCtrl,"CFLRestart",cflv[0]);
00077     RegisterAt(LocCtrl,"CFLControl",cflv[1]);
00078     RegisterAt(LocCtrl,"Outputs",Outputs);
00079   }
00080   
00081 public: 
00082   double LastTime;
00083   int VariableTimeStepping;
00084   int Outputs;
00085   double dtv[2], cflv[2];
00086   ControlDevice LocCtrl;     
00087 };
00088 
00089 
00099 template <class VectorType, class FixupType, class FlagType>
00100 class SolverControl(DIM) : public SolverControlInterface(DIM)<VectorType> {
00101   typedef Solver(DIM)<VectorType,FixupType,FlagType> solver_type;
00102   typedef ExactSolution(DIM)<VectorType> exact_solution_type;
00103   typedef typename VectorType::InternalDataType DataType;
00104 
00105   typedef GridFunction(DIM)<VectorType> vec_grid_fct_type;  
00106   typedef GridFunction(DIM)<DataType>   grid_fct_type;   
00107 
00108 public:
00109   SolverControl(DIM)(solver_type& solver) : 
00110     StepControls(1),  OutEvery(0), CheckPointEvery(0), 
00111     RedistributeEvery(0), InitOutput(1), Restart(0), Distribution(DAGHCompositeDistribution), 
00112     MaxLev(1), GuCFactor(2), CutOffs(0), _Solver(solver), ENorm(0) {    
00113     for (register int d=0; d<DIM; d++) {
00114       shape[d] = 2;
00115       geom[2*d] = 0.0;
00116       geom[2*d+1] = 1.0;
00117       periodic[d] = DAGHFalse;
00118     }    
00119     for (register int i=0; i<DAGHMaxLevels; i++)
00120       RefineFactor[i] = DAGHDefaultRefineFactor;
00121     CompName = (string*) 0;
00122     GH = (GridHierarchy *) 0;
00123     _ExactSolution = (exact_solution_type *) 0;
00124     CheckpointName = "checkpt";
00125   }
00126 
00127   ~SolverControl(DIM)() {
00128     if (CompName) delete [] CompName;
00129     if (GH) delete GH;
00130   }
00131  
00132   virtual void register_at(ControlDevice& Ctrl) {}
00133   virtual void register_at(ControlDevice& Ctrl, const string& prefix) {
00134     LocCtrl = Ctrl.getSubDevice(prefix+"SolverControl");
00135     char VariableName[32];
00136     for (int d=0; d<DIM; d++) {
00137       sprintf(VariableName,"Cells(%d)",d+1);
00138       RegisterAt(LocCtrl,VariableName,shape[d]);
00139       sprintf(VariableName,"GeomMin(%d)",d+1);
00140       RegisterAt(LocCtrl,VariableName,geom[2*d]);
00141       sprintf(VariableName,"GeomMax(%d)",d+1);
00142       RegisterAt(LocCtrl,VariableName,geom[2*d+1]);     
00143       sprintf(VariableName,"PeriodicBoundary(%d)",d+1);
00144       RegisterAt(LocCtrl,VariableName,periodic[d]);     
00145     }
00146     for (int nc=0; nc<MaxCutOffRegions; nc++) {
00147       for (int d=0; d<DIM; d++) 
00148         for (int j=0; j<2; j++) {
00149           sprintf(VariableName,"CCells(%d,%d,%d)",nc+1,j+1,d+1);
00150           RegisterAt(LocCtrl,VariableName,cuts[nc*2*DIM+2*d+j]); 
00151           cuts[nc*2*DIM+2*d+j] = -100000;
00152         }
00153     }
00154     RegisterAt(LocCtrl,"CutOffs",CutOffs); 
00155     for (register int i=0; i<DAGHMaxLevels; i++) {
00156       sprintf(VariableName,"RefineFactor(%d)",i+1);
00157       RegisterAt(LocCtrl,VariableName,RefineFactor[i]);
00158     }
00159     for (int cs=0; cs<MaxStepControls; cs++) {
00160       char text[5];
00161       sprintf(text,"%d-",cs+1);
00162       Step[cs].register_at(LocCtrl,string(text));
00163     }
00164     RegisterAt(LocCtrl,"StepControls",StepControls); 
00165     RegisterAt(LocCtrl,"CheckEvery",CheckPointEvery);
00166     RegisterAt(LocCtrl,"RedistributeEvery",RedistributeEvery);
00167     RegisterAt(LocCtrl,"Restart",Restart);  
00168     RegisterAt(LocCtrl,"InitOutput",InitOutput);      
00169     RegisterAt(LocCtrl,"MaxLevels",MaxLev);
00170     RegisterAt(LocCtrl,"Distribution",Distribution);
00171     RegisterAt(LocCtrl,"CheckpointName",CheckpointName);    
00172     RegisterAt(LocCtrl,"ErrorNorm",ENorm);    
00173     RegisterAt(LocCtrl,"GuCFactor",GuCFactor);    
00174     
00175     if (CompName) delete [] CompName;
00176     CompName = new string[NOutputs()];
00177     for (int comp=0; comp<NOutputs(); comp++) {
00178       char ComponentName[16];
00179       sprintf(ComponentName,"-");
00180       CompName[comp] = ComponentName;
00181       sprintf(ComponentName,"OutputName(%d)",comp+1);
00182       RegisterAt(LocCtrl,ComponentName,CompName[comp]);
00183     }    
00184   }
00185 
00186   virtual void init() {
00187     DAGHIOInit();   
00188     Solver_().init();
00189   }
00190 
00191   virtual void update() {    
00192     Solver_().update();
00193   }
00194 
00195   virtual void finish() {
00196     DAGHIO_HDF_NCSA_Finalize();   
00197 #ifndef DAGH_NO_MPI
00198     cout << "Finalizing MPI\n" << flush;
00199     DAGHMPI_Finalize();
00200 #endif  
00201     Solver_().finish();
00202     if (CompName) {
00203       delete [] CompName;
00204       CompName = (string*) 0;
00205     }
00206     if (GH) {
00207       delete GH;
00208       GH = (GridHierarchy *) 0;
00209     }
00210   } 
00211 
00212   //---------------------------------------------------------------------- 
00213   // Output in conservative variables   
00214   //---------------------------------------------------------------------- 
00215   virtual void WriteOut(GridHierarchy& GH, vec_grid_fct_type& u,
00216                         grid_fct_type& IOfunc) {
00217     int Time;       
00218     int me = MY_PROC(GH); 
00219     
00220     Time = CurrentTime(GH,0); 
00221     Solver_().SetLastOutputTime(Time);
00222     char ioname[DAGHBktGFNameWidth+16];
00223     for (int i=0; i<NEquations(); i++) { 
00224       if (CompName[i].c_str()[0] == '-') 
00225         continue;
00226       sprintf(ioname,"%s_%d",CompName[i].c_str(),Time); 
00227       if (me == VizServer) 
00228         cout << " *** Writing " << ioname << endl;   
00229       for (int lev=0; lev<=FineLevel(GH); lev++) {
00230         Time = CurrentTime(GH,lev); 
00231         forall (u,Time,lev,c)           
00232           equals_from(IOfunc(Time,lev,c), u(Time,lev,c), i);
00233         end_forall      
00234         SetPhysicalTime(IOfunc,Time,lev,GetPhysicalTime(u,Time,lev));
00235         Write(IOfunc,Time,lev,DAGH_Double,ioname);      
00236       }
00237       // This barrier ensures that SDSclose() is called only once!
00238 #ifndef DAGH_NO_MPI
00239       MPI_Barrier(comm_service::comm()); 
00240 #endif
00241     } 
00242   }   
00243 
00244   //---------------------------------------------------------------------- 
00245   // Calculate Error Norm for conservative variables 
00246   //---------------------------------------------------------------------- 
00247   virtual void ErrorNorm(GridHierarchy& GH, vec_grid_fct_type& u, double t,
00248                          grid_fct_type& work) {
00249     if (ENorm<1 || ENorm>3) 
00250       return;
00251 
00252     int me = MY_PROC(GH);
00253     if (me == VizServer) 
00254       cout << endl << "                         Errors:";
00255     for (int lev=0; lev<=FineLevel(GH); lev++) {
00256       int Time = CurrentTime(GH,lev); 
00257       int TStep = TimeStep(u,lev);
00258       DCoords dx(DIM);
00259       for(int d=0; d<DIM; d++)
00260         dx(d) = DeltaX(u,d,lev);
00261       
00262       forall (u,Time+TStep,lev,c) 
00263         TheExactSolution().SetGrid(u(Time+TStep,lev,c), t, dx);
00264         u(Time+TStep,lev,c).minus(u(Time,lev,c));
00265       end_forall      
00266 
00267       double Error = Norm(u,Time+TStep,lev,ENorm);
00268       if (me == VizServer) 
00269         cout << "   L(" << lev << ")=" << Error;
00270     }
00271   }  
00272 
00273   virtual void Solve() {
00274     START_WATCH_WHOLE
00275     //---------------------------------------------------------------
00276     // Setup grid hierarchy
00277     //---------------------------------------------------------------
00278     GH = new GridHierarchy(DIM,DAGHCellCentered,MaxLev);
00279     //---------------------------------------------------------------
00280     // Processor information
00281     //---------------------------------------------------------------
00282     int me = MY_PROC(*GH);    
00283     //---------------------------------------------------------------
00284     // Dimensioning of the grid. Avoid negative indices!
00285     //---------------------------------------------------------------
00286     // for (int d=0; d<DIM; d++)
00287     //   offset[d] = GH.daghshadowfactor()*NGhosts();
00288  
00289     SetBaseGrid(*GH,geom,shape,CutOffs,cuts);   
00290     for (register int l=0; l<MaxLev-1; l++) {
00291       if (RefineFactor[l]<DAGHDefaultRefineFactor)
00292         RefineFactor[l] = DAGHDefaultRefineFactor;
00293       if (Solver_().ErrorEstimation() && RefineFactor[l]%DAGHShadowFactor) {
00294         cout << "   RefineFactor(" << l+1 << ") is not dividable by Shadowfactor()." 
00295              << " Can't use error estimation.\n" << flush;
00296         Solver_().Flagging_().SetErrorEstimation(false);
00297       } 
00298       SetRefineFactor(*GH,l,RefineFactor[l]); 
00299     }
00300     SetBoundaryWidth(*GH,NGhosts());
00301     SetBoundaryType(*GH,DAGHBoundaryUserDef); 
00302     if (Solver_().ErrorEstimation()) 
00303       GuCFactor = Max(GuCFactor,DAGHShadowFactor);
00304     if (GuCFactor) SetGuCFactor(*GH,GuCFactor);
00305     for (int d=0; d<DIM; d++)
00306       if (periodic[d]) SetPeriodicBoundaries(*GH,d,DAGHTrue);
00307     SetDistributionType(*GH,Distribution);
00308     
00309     DAGHIOType(*GH, DAGHIO_HDF_NCSA);   
00310     ComposeHierarchy(*GH);     
00311 
00312     Solver_().SetGridHierarchy(GH);
00313     if (_ExactSolution) 
00314       TheExactSolution().SetGridHierarchy(GH);
00315 
00316     BEGIN_COMPUTE               
00317       //---------------------------------------------------------------
00318       //  Declare grid functions    
00319       //---------------------------------------------------------------
00320       AllocError::SetTexts("main()","allocation of GridFunctions");
00321 
00322       int t_sten = 1;   
00323       int s_sten = NGhosts();       
00324       char GFName[DAGHBktGFNameWidth], GFNamesh[DAGHBktGFNameWidth];
00325       sprintf(GFName,"u");
00326       vec_grid_fct_type u(GFName, t_sten, s_sten, *GH, DAGHCommSimple);
00327       SetTimeAlias(u, -1, 1);
00328       
00329       vec_grid_fct_type* ush = (vec_grid_fct_type*) 0;
00330       if (Solver_().ErrorEstimation()) {
00331         sprintf(GFNamesh,"ush");        
00332         ush = new vec_grid_fct_type(GFNamesh, t_sten, s_sten,  
00333                                     DAGHShadowFactor, *GH, DAGHCommSimple);
00334         SetTimeAlias(*ush, -1, 1);
00335       }
00336       
00337       // One work array with two time steps and Shadow-hierarchy is necessary - 
00338       // Could be reused in AMRClawFlagging(DIM) 
00339       char IOName[DAGHBktGFNameWidth], IONamesh[DAGHBktGFNameWidth];
00340       sprintf(IOName,"IOPlaceholder"); 
00341       grid_fct_type IOfunc(IOName, t_sten, s_sten, *GH, DAGHCommSimple);
00342       SetCheckpointFlag(IOfunc, DAGHFalse);
00343       IOfunc.GF_SetMaxRecomposeLevel(DAGHNull);
00344       SetTimeAlias(IOfunc, -1, 1);
00345  
00346       grid_fct_type* IOfuncsh = (grid_fct_type*) 0;
00347       if (Solver_().ErrorEstimation()) {
00348         sprintf(IONamesh,"IOPlaceholdersh"); 
00349         IOfuncsh = new grid_fct_type(IONamesh, t_sten, s_sten,  
00350                                      DAGHShadowFactor, *GH, DAGHCommSimple);
00351         SetCheckpointFlag(*IOfuncsh, DAGHFalse);
00352         IOfuncsh->GF_SetMaxRecomposeLevel(DAGHNull);
00353         SetTimeAlias(*IOfuncsh, -1, 1);
00354       }
00355  
00356       Solver_().SetGridFunctions(&u, ush, &IOfunc, IOfuncsh);
00357       Solver_().SetupData();      
00358       Solver_().SetBoundaryConditions();
00359                   
00360       //---------------------------------------------------------------
00361       // Begin Evolution        
00362       //---------------------------------------------------------------    
00363       int CStepControl = 0;    
00364       double Current_Time, Next_Time;
00365       if (!Restart) {
00366         if (me == VizServer)    
00367           cout << " --- Initializing --- " << flush;
00368         Solver_().Initialize(Step[0].dtv[0]);     
00369         if (me == VizServer)       
00370           cout << "   Levels = " << FineLevel(*GH)+1;
00371       } 
00372       else {  
00373         if (me == VizServer)     
00374           cout << " --- Restarting from " << CheckpointName.c_str() << " --- " << flush;
00375         Solver_().Restart(CheckpointName.c_str());
00376         if (me == VizServer)    
00377           cout << "   Last Iteration: " << StepsTaken(*GH,0);
00378       }  
00379       GH->DAGH_GetTimeSpecs(Current_Time, Next_Time);
00380       if (_ExactSolution) 
00381         ErrorNorm(*GH, u, Current_Time, IOfunc);
00382       if (me == VizServer)    
00383         cout << endl << flush;
00384       while (Current_Time>=Step[CStepControl].LastTime &&
00385              CStepControl<StepControls)
00386         CStepControl++;  
00387         
00388       int FirstStep = StepsTaken(*GH,0);
00389       double tout = (CStepControl>0 ? Step[CStepControl-1].LastTime : 0.0);
00390       double dtout = Step[CStepControl].LastTime - 
00391         (CStepControl>0 ? Step[CStepControl-1].LastTime : 0.0);
00392       if (Step[CStepControl].Outputs > 0)
00393           dtout /= Step[CStepControl].Outputs;
00394       while (Current_Time+DBL_EPSILON*FACTOR > tout+dtout)
00395         tout += dtout;
00396 
00397       while (Step[StepControls-1].LastTime-tout > DBL_EPSILON*FACTOR) {   
00398         if (Step[CStepControl].Outputs>0 && InitOutput) {
00399           START_WATCH 
00400             WriteOut(*GH, u, IOfunc);     
00401           END_WATCH_OUPUT  
00402         }
00403         InitOutput = 1;
00404 
00405         while (tout+dtout-Current_Time > DBL_EPSILON*FACTOR) {   
00406           if (CheckPointEvery > 0 && StepsTaken(*GH,0) > FirstStep) 
00407             if (StepsTaken(*GH,0) % CheckPointEvery == 0) {
00408               if (me == VizServer) 
00409                 cout << " *** Writing Checkpoint-File" << endl << flush;   
00410               START_WATCH 
00411 #ifndef DAGH_NO_MPI
00412                 MPI_Barrier(comm_service::comm()); 
00413 #endif
00414                 DAGHIO_HDF_NCSA_Finalize();           
00415                 Checkpoint(*GH, CheckpointName.c_str());
00416 #ifndef DAGH_NO_MPI
00417                 MPI_Barrier(comm_service::comm()); 
00418 #endif
00419               END_WATCH_OUPUT  
00420             }       
00421           
00422           if (Next_Time > tout+dtout)
00423             if (tout+dtout-Current_Time > DBL_EPSILON*FACTOR)
00424               SetTimeSpecs(*GH, Current_Time, tout+dtout, 2);
00425             else
00426               break;
00427           
00428           if (me == VizServer) { 
00429             cout << " --- Iteration: " << StepsTaken(*GH,0)+1 << " ---";
00430             cout << "   t = " << Current_Time << flush;  
00431           } 
00432     
00433           int Rejections;
00434           double cfl_max = Solver_().Tick(Step[CStepControl].VariableTimeStepping, 
00435                                           Step[CStepControl].dtv, Step[CStepControl].cflv,
00436                                           Rejections,RedistributeEvery);
00437  
00438           double Help_Time;
00439           GH->DAGH_GetTimeSpecs(Next_Time, Help_Time);
00440           if (me == VizServer) { 
00441             cout << "  dt=" << Next_Time-Current_Time 
00442                  << "  Levels=" << FineLevel(*GH)+1 
00443                  << "  max. CFL=" << cfl_max;
00444           }
00445           Current_Time = Next_Time;
00446           Next_Time = Help_Time;
00447           
00448           if (_ExactSolution) 
00449             ErrorNorm(*GH, u, Current_Time, IOfunc);
00450           if (me == VizServer) {
00451             if (Rejections)
00452               cout << "  Restarted: " << Rejections << "x";
00453             cout << endl << flush;
00454           }
00455         }  
00456         tout += dtout;
00457         if (tout >= Step[CStepControl].LastTime) {
00458           CStepControl++;     
00459           if (CStepControl<StepControls) {   
00460             dtout = Step[CStepControl].LastTime - Step[CStepControl-1].LastTime;
00461             if (Step[CStepControl].Outputs > 0)
00462               dtout /= Step[CStepControl].Outputs;
00463           }
00464         }
00465       }
00466       START_WATCH 
00467         WriteOut(*GH, u, IOfunc);
00468       END_WATCH_OUPUT   
00469 
00470       if (CheckPointEvery > 0) {
00471         if (me == VizServer) 
00472           cout << " *** Writing Checkpoint-File" << endl << flush;   
00473         START_WATCH 
00474 #ifndef DAGH_NO_MPI
00475           MPI_Barrier(comm_service::comm()); 
00476 #endif
00477           DAGHIO_HDF_NCSA_Finalize();                 
00478           Checkpoint(*GH, CheckpointName.c_str());
00479 #ifndef DAGH_NO_MPI
00480           MPI_Barrier(comm_service::comm()); 
00481 #endif
00482         END_WATCH_OUPUT         
00483       }
00484   
00485       if (ush) delete ush;
00486       if (IOfuncsh) delete IOfuncsh;     
00487 
00488       START_WATCH 
00489 #ifndef DAGH_NO_MPI
00490         double start_finalize = MPI_Wtime();
00491         while (MPI_Wtime()-start_finalize < 10.0) {}
00492         MPI_Barrier(comm_service::comm()); 
00493 #endif
00494         DAGHIOEnd(*GH);  
00495       END_WATCH_OUPUT   
00496 
00497     END_COMPUTE
00498     END_WATCH_WHOLE
00499 
00500     if (me == VizServer)    
00501       Timing::print(cout);
00502   } 
00503   inline void SetExactSolution(exact_solution_type* exact) { _ExactSolution = exact; }
00504   inline exact_solution_type& TheExactSolution() { return *_ExactSolution; }
00505   inline const exact_solution_type& TheExactSolution() const { return *_ExactSolution; }
00506 
00507   inline solver_type& Solver_() { return _Solver; }
00508   inline const solver_type& Solver_() const { return _Solver; }
00509   inline const int& NEquations() const { return _Solver.NEquations(); }
00510   inline const int& NGhosts() const { return _Solver.NGhosts(); }
00511   virtual int NOutputs() const { return NEquations(); }
00512 
00513 private:
00514   int StepControls;             
00515   StepControl Step[MaxStepControls];
00516   int OutEvery;         
00517   int CheckPointEvery;
00518   int RedistributeEvery;
00519   int InitOutput;
00520   int Restart;          
00521   int Distribution;  
00522   int MaxLev;   
00523   int GuCFactor;
00524   int RefineFactor[DAGHMaxLevels];
00525   int periodic[DIM];
00526   int shape[DIM];  
00527   double geom[2*DIM];
00528   int CutOffs;
00529   int cuts[2*DIM*MaxCutOffRegions];
00530   char chkptdata[DAGHChkPtTagNameSize];
00531   string CheckpointName;
00532   GridHierarchy* GH;
00533 
00534 protected:
00535   string* CompName;
00536   solver_type& _Solver; 
00537   exact_solution_type* _ExactSolution;
00538   int ENorm;
00539   ControlDevice LocCtrl;
00540 };
00541 
00542 
00543 #endif
00544
SolverControl.h

last update: 06/01/04
