PISM, A Parallel Ice Sheet Model: src/base/iceModel.hh Source File

PISM, A Parallel Ice Sheet Model stable 0.4.1779
00001 // Copyright (C) 2004-2011 Jed Brown, Ed Bueler and Constantine Khroulev
00002 //
00003 // This file is part of PISM.
00004 //
00005 // PISM is free software; you can redistribute it and/or modify it under the
00006 // terms of the GNU General Public License as published by the Free Software
00007 // Foundation; either version 2 of the License, or (at your option) any later
00008 // version.
00009 //
00010 // PISM is distributed in the hope that it will be useful, but WITHOUT ANY
00011 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00012 // FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
00013 // details.
00014 //
00015 // You should have received a copy of the GNU General Public License
00016 // along with PISM; if not, write to the Free Software
00017 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
00018 
00019 #ifndef __iceModel_hh
00020 #define __iceModel_hh
00021 
00023 
00036 #include <signal.h>
00037 #include <gsl/gsl_rng.h>
00038 #include <petscsnes.h>
00039 #include <petsctime.h>          // PetscGetTime()
00040 
00041 #include "flowlaw_factory.hh"   // IceFlowLawFactory and friends
00042 #include "basal_resistance.hh"  // IceBasalResistancePlasticLaw
00043 #include "PISMYieldStress.hh"
00044 #include "pism_const.hh"
00045 #include "enthalpyConverter.hh"
00046 #include "grid.hh"
00047 #include "iceModelVec.hh"
00048 #include "NCVariable.hh"
00049 #include "PISMBedSmoother.hh"
00050 #include "PISMVars.hh"
00051 #include "Timeseries.hh"
00052 #include "PISMStressBalance.hh"
00053 #include "bedrockThermalUnit.hh"
00054 
00055 #include "PISMBedDef.hh"
00056 #include "PISMOcean.hh"
00057 #include "PISMSurface.hh"
00058 
00059 // use namespace std BUT remove trivial namespace browser from doxygen-erated HTML source browser
00061 using namespace std;
00063 
00064 
00066 class IceModel {
00067   // The following classes implement various diagnostic computations.
00068   friend class IceModel_hardav;
00069   friend class IceModel_bwp;
00070   friend class IceModel_cts;
00071   friend class IceModel_dhdt;
00072   friend class IceModel_temp;
00073   friend class IceModel_temp_pa;
00074   friend class IceModel_temppabase;
00075   friend class IceModel_enthalpybase;
00076   friend class IceModel_enthalpysurf;
00077   friend class IceModel_tempbase;
00078   friend class IceModel_tempsurf;
00079   friend class IceModel_liqfrac;
00080   friend class IceModel_tempicethk;
00081   friend class IceModel_tempicethk_basal;
00082   friend class IceModel_new_mask;
00083 public:
00084   // see iceModel.cc for implementation of constructor and destructor:
00085   IceModel(IceGrid &g, NCConfigVariable &config, NCConfigVariable &overrides);
00086   virtual ~IceModel(); // must be virtual merely because some members are virtual
00087 
00088   // see iMinit.cc
00089   virtual PetscErrorCode grid_setup();
00090 
00091   virtual PetscErrorCode allocate_submodels();
00092   virtual PetscErrorCode allocate_flowlaw();
00093   virtual PetscErrorCode allocate_enthalpy_converter();
00094   virtual PetscErrorCode allocate_basal_resistance_law();
00095   virtual PetscErrorCode allocate_stressbalance();
00096   virtual PetscErrorCode allocate_bed_deformation();
00097   virtual PetscErrorCode allocate_bedrock_thermal_unit();
00098   virtual PetscErrorCode allocate_basal_yield_stress();
00099 
00100   virtual PetscErrorCode init_couplers();
00101   virtual PetscErrorCode set_grid_from_options();
00102   virtual PetscErrorCode set_grid_defaults();
00103   virtual PetscErrorCode model_state_setup();
00104   virtual PetscErrorCode set_vars_from_options();
00105   virtual PetscErrorCode allocate_internal_objects();
00106   virtual PetscErrorCode misc_setup();
00107   virtual PetscErrorCode init_diagnostics();
00108   virtual PetscErrorCode init_ocean_kill();
00109 
00110   // see iceModel.cc
00111   PetscErrorCode init();
00112   virtual PetscErrorCode run();
00113   virtual PetscErrorCode step(bool do_mass_continuity, 
00114                               bool do_energy,
00115                               bool do_age,
00116                               bool do_skip,
00117                               bool use_ssa_when_grounded);
00118   virtual PetscErrorCode setExecName(const char *my_executable_short_name);
00119   virtual IceFlowLawFactory &getIceFlowLawFactory() { return iceFactory; }
00120   virtual IceFlowLaw *getIceFlowLaw() {return ice;}
00121   virtual void reset_counters();
00122 
00123   // see iMbootstrap.cc 
00124   virtual PetscErrorCode bootstrapFromFile(const char *fname);
00125   virtual PetscErrorCode bootstrap_2d(const char *fname);
00126   virtual PetscErrorCode bootstrap_3d();
00127   virtual PetscErrorCode putTempAtDepth();
00128 
00129   // see iMoptions.cc
00130   virtual PetscErrorCode setFromOptions();
00131   virtual PetscErrorCode set_output_size(string option, string description,
00132                                          string default_value, set<string> &result);
00133 
00134   // see iMutil.cc
00135   virtual void attach_surface_model(PISMSurfaceModel *surf);
00136   virtual void attach_ocean_model(PISMOceanModel *ocean);
00137   virtual PetscErrorCode additionalAtStartTimestep();
00138   virtual PetscErrorCode additionalAtEndTimestep();
00139   virtual PetscErrorCode compute_cell_areas(); // is an initialization step; should go there
00140 
00141   // see iMIO.cc
00142   virtual PetscErrorCode initFromFile(const char *);
00143   virtual PetscErrorCode writeFiles(const char* default_filename);
00144   virtual PetscErrorCode write_model_state(const char *filename);
00145   virtual PetscErrorCode write_metadata(const char *filename);
00146   virtual PetscErrorCode write_variables(const char* filename, set<string> vars,
00147                                          nc_type nctype);
00148 protected:
00149 
00150   IceGrid               &grid;
00151 
00152   NCConfigVariable      mapping, 
00153     &config,                     
00154     &overrides;                  
00155   NCGlobalAttributes    global_attributes;
00156 
00157   IceFlowLawFactory     iceFactory;
00158   IceFlowLaw            *ice;
00159 
00160   PISMYieldStress *basal_yield_stress;
00161   IceBasalResistancePlasticLaw *basal;
00162 
00163   EnthalpyConverter *EC;
00164   PISMBedThermalUnit *btu;
00165 
00166   PISMSurfaceModel *surface;
00167   PISMOceanModel   *ocean;
00168   PISMBedDef       *beddef;
00169 
00172   PISMVars variables;
00173 
00174   // state variables and some diagnostics/internals
00175   IceModelVec2S
00176         vh,             
00177         vH,             
00178         vdHdt,          
00179         vtauc,          
00180         vHmelt,         
00181         vbmr,       
00182         vLongitude,     
00183         vLatitude,      
00184         vbed,           
00185         vuplift,        
00186         vGhf,           
00187         bedtoptemp,     
00188 
00189 
00190         vHref,          
00191         vHresidual,     
00192         vPrinStrain1,   
00193         vPrinStrain2,   
00194 
00195     acab,               
00196     artm,               
00197     liqfrac_surface,    
00198     shelfbtemp,         
00199     shelfbmassflux,     
00200         cell_area;              
00201 
00202         
00203  
00204   IceModelVec2Int vMask, 
00205 
00206     ocean_kill_mask,     
00207     vIcebergMask, 
00208 
00209         vBCMask; 
00210  
00211   IceModelVec2V vBCvel; 
00212 
00213 
00214   IceModelVec3
00215         T3,             
00216         Enth3,          
00217         tau3;           
00218 
00219   // parameters
00220   PetscReal   dt,     
00221               t_years_TempAge,  
00222               dt_years_TempAge,  
00223               maxdt_temporary, dt_force,
00224               CFLviolcount,    
00225               dt_from_diffus, dt_from_cfl, CFLmaxdt, CFLmaxdt2D,
00226               gDmax,            // global max of the diffusivity
00227               gmaxu, gmaxv, gmaxw,  // global maximums on 3D grid of abs value of vel components
00228               gdHdtav,  
00229     total_sub_shelf_ice_flux,
00230     total_basal_ice_flux,
00231     total_surface_ice_flux,
00232     nonneg_rule_flux,
00233     ocean_kill_flux,
00234     float_kill_flux,
00235     dvoldt;  
00236   PetscInt    skipCountDown;
00237 
00238   // physical parameters used frequently enough to make looking up via
00239   // config.get() a hassle:
00240   PetscScalar standard_gravity;
00241   // Initialized in the IceModel constructor from the configuration file;
00242   // SHOULD NOT be hard-wired.
00243 
00244   // flags
00245   PetscTruth  updateHmelt, shelvesDragToo, allowAboveMelting;
00246   PetscTruth  repeatRedist, putOnTop;
00247   char        adaptReasonFlag;
00248 
00249   string      stdout_flags, stdout_ssa;
00250 
00251   string executable_short_name;
00252   
00253 protected:
00254   // see iceModel.cc
00255   virtual PetscErrorCode createVecs();
00256   virtual PetscErrorCode deallocate_internal_objects();
00257 
00258   // see iMadaptive.cc
00259   virtual PetscErrorCode computeMax3DVelocities();
00260   virtual PetscErrorCode computeMax2DSlidingSpeed();
00261   virtual PetscErrorCode adaptTimeStepDiffusivity();
00262   virtual PetscErrorCode determineTimeStep(const bool doTemperatureCFL);
00263   virtual PetscErrorCode countCFLViolations(PetscScalar* CFLviol);
00264 
00265   // see iMage.cc
00266   virtual PetscErrorCode ageStep();
00267 
00268   // see iMbeddef.cc
00269   PetscScalar last_bed_def_update;
00270   virtual PetscErrorCode bed_def_step(bool &bed_changed);
00271 
00272   // see iMcalving.cc
00273   virtual PetscErrorCode eigenCalving();
00274   virtual PetscErrorCode calvingAtThickness();
00275 
00276   // see iMdefaults.cc
00277   PetscErrorCode setDefaults();
00278 
00279   // see iMenergy.cc
00280   virtual PetscErrorCode energyStep();
00281   virtual PetscErrorCode get_bed_top_temp(IceModelVec2S &result);
00282   virtual bool checkThinNeigh(
00283        PetscScalar E, PetscScalar NE, PetscScalar N, PetscScalar NW, 
00284        PetscScalar W, PetscScalar SW, PetscScalar S, PetscScalar SE);
00285 
00286   // see iMenthalpy.cc
00287   virtual PetscErrorCode compute_enthalpy_cold(IceModelVec3 &temperature, IceModelVec3 &result);
00288   virtual PetscErrorCode compute_enthalpy(IceModelVec3 &temperature, IceModelVec3 &liquid_water_fraction,
00289                                           IceModelVec3 &result);
00290   virtual PetscErrorCode compute_liquid_water_fraction(IceModelVec3 &enthalpy, IceModelVec3 &result);
00291 
00292   virtual PetscErrorCode setCTSFromEnthalpy(IceModelVec3 &useForCTS);
00293 
00294   virtual PetscErrorCode getEnthalpyCTSColumn(PetscScalar p_air, 
00295                                               PetscScalar thk,   
00296                                               PetscInt ks,       
00297                                               const PetscScalar *Enth, 
00298                                               const PetscScalar *w,    
00299                                               PetscScalar *lambda,     
00300                                               PetscScalar **Enth_s
00301                                               );   
00302   virtual PetscErrorCode enthalpyAndDrainageStep(
00303                 PetscScalar* vertSacrCount, PetscScalar* liquifiedVol,
00304                 PetscScalar* bulgeCount);
00305 
00306   // see iMgeometry.cc
00307   virtual PetscErrorCode updateSurfaceElevationAndMask();
00308   virtual PetscErrorCode update_mask();
00309   virtual PetscErrorCode update_surface_elevation();
00310   virtual PetscErrorCode massContExplicitStep();
00311 
00312   // see iMicebergs.cc
00313   virtual PetscErrorCode killIceBergs();           // call this one to do proper sequence
00314   virtual PetscErrorCode findIceBergCandidates();
00315   virtual PetscErrorCode identifyNotAnIceBerg();
00316   virtual PetscErrorCode killIdentifiedIceBergs();
00317   virtual PetscErrorCode killEasyIceBergs();       // FIXME: do we want this one to happen even if eigencalving does not happen?  should we be calling this one before any time that principal values need to be computed?
00318 
00319   // see iMIO.cc
00320   virtual PetscErrorCode set_time_from_options();
00321   virtual PetscErrorCode dumpToFile(const char *filename);
00322   virtual PetscErrorCode regrid(int dimensions);
00323   virtual PetscErrorCode regrid_variables(string filename, set<string> regrid_vars, int ndims);
00324 
00325   // see iMpartgrid.cc
00326   PetscErrorCode cell_interface_velocities(bool do_part_grid,
00327                                            int i, int j,
00328                                            planeStar<PetscScalar> &vel_output);
00329   PetscReal get_average_thickness(bool do_redist, planeStar<int> M,
00330                                   planeStar<PetscScalar> H);
00331   virtual PetscErrorCode redistResiduals();
00332   virtual PetscErrorCode calculateRedistResiduals();
00333 
00334   // see iMreport.cc
00335   virtual PetscErrorCode volumeArea(
00336                        PetscScalar& gvolume,PetscScalar& garea,
00337                        PetscScalar& gvolSIA, PetscScalar& gvolstream, 
00338                        PetscScalar& gvolshelf);
00339   virtual PetscErrorCode energyStats(
00340                        PetscScalar iarea,
00341                        PetscScalar &gmeltfrac, PetscScalar &gtemp0);
00342   virtual PetscErrorCode ageStats(PetscScalar ivol, PetscScalar &gorigfrac);
00343   virtual PetscErrorCode summary(bool tempAndAge);
00344   virtual PetscErrorCode summaryPrintLine(
00345               PetscTruth printPrototype, bool tempAndAge,
00346               PetscScalar year, PetscScalar delta_t, 
00347               PetscScalar volume, PetscScalar area,
00348               PetscScalar meltfrac, PetscScalar H0, PetscScalar T0);
00349 
00350   // see iMreport.cc;  methods for computing diagnostic quantities:
00351   // scalar:
00352   virtual PetscErrorCode ice_mass_bookkeeping();
00353   virtual PetscErrorCode compute_ice_volume(PetscScalar &result);
00354   virtual PetscErrorCode compute_ice_volume_temperate(PetscScalar &result);
00355   virtual PetscErrorCode compute_ice_volume_cold(PetscScalar &result);
00356   virtual PetscErrorCode compute_ice_area(PetscScalar &result);
00357   virtual PetscErrorCode compute_ice_area_temperate(PetscScalar &result);
00358   virtual PetscErrorCode compute_ice_area_cold(PetscScalar &result);
00359   virtual PetscErrorCode compute_ice_area_grounded(PetscScalar &result);
00360   virtual PetscErrorCode compute_ice_area_floating(PetscScalar &result);
00361   virtual PetscErrorCode compute_ice_enthalpy(PetscScalar &result);
00362   virtual PetscErrorCode compute_by_name(string name, PetscScalar &result);
00363   
00364   // see iMtemp.cc
00365   virtual PetscErrorCode excessToFromBasalMeltLayer(
00366                       const PetscScalar rho_c, const PetscScalar z, const PetscScalar dz,
00367                       PetscScalar *Texcess, PetscScalar *Hmelt);
00368   virtual PetscErrorCode temperatureStep(PetscScalar* vertSacrCount, PetscScalar* bulgeCount);
00369 
00370   // see iMutil.cc
00371   virtual int            endOfTimeStepHook();
00372   virtual PetscErrorCode stampHistoryCommand();
00373   virtual PetscErrorCode stampHistoryEnd();
00374   virtual PetscErrorCode stampHistory(string);
00375   virtual PetscErrorCode check_maximum_thickness();
00376   virtual PetscErrorCode check_maximum_thickness_hook(const int old_Mz);
00377   virtual bool           issounding(const PetscInt i, const PetscInt j);
00378 
00379 protected:
00380   // working space (a convenience)
00381   static const PetscInt nWork2d=2;
00382   IceModelVec2S vWork2d[nWork2d];
00383   IceModelVec2V vWork2dV;
00384 
00385   // 3D working space
00386   IceModelVec3 vWork3d;
00387 
00388   PISMStressBalance *stress_balance;
00389 
00390   map<string,PISMDiagnostic*> diagnostics;
00391 
00392   // Set of variables to put in the output file:
00393   set<string> output_vars;
00394 
00395   // This is related to the snapshot saving feature
00396   string snapshots_filename;
00397   bool save_snapshots, snapshots_file_is_ready, split_snapshots;
00398   vector<double> snapshot_times;
00399   set<string> snapshot_vars;
00400   unsigned int current_snapshot;
00401   PetscErrorCode init_snapshots();
00402   PetscErrorCode write_snapshot();
00403 
00404   // scalar time-series
00405   bool save_ts;                 
00406   string ts_filename;           
00407   vector<double> ts_times;      
00408   unsigned int current_ts;      
00409   set<string> ts_vars;          
00410   vector<DiagnosticTimeseries*> timeseries;
00411   PetscErrorCode init_timeseries();
00412   PetscErrorCode create_timeseries();
00413   PetscErrorCode flush_timeseries();
00414   PetscErrorCode write_timeseries();
00415   PetscErrorCode ts_max_timestep(double t_years, double& dt_years);
00416 
00417   // spatially-varying time-series
00418   bool save_extra, extra_file_is_ready, split_extra;
00419   string extra_filename;
00420   vector<double> extra_times;
00421   unsigned int current_extra;
00422   set<string> extra_vars;
00423   PetscErrorCode init_extras();
00424   PetscErrorCode write_extras();
00425   PetscErrorCode extras_max_timestep(double t_years, double& dt_years);
00426 
00427   // automatic backups
00428   double backup_interval;
00429   string backup_filename;
00430   PetscReal last_backup_time;
00431   set<string> backup_vars;
00432   PetscErrorCode init_backups();
00433   PetscErrorCode write_backup();
00434 
00435   // diagnostic viewers; see iMviewers.cc
00436   virtual PetscErrorCode init_viewers();
00437   virtual PetscErrorCode update_viewers();
00438   set<string> map_viewers, slice_viewers, sounding_viewers;
00439   PetscInt     id, jd;       // sounding indices
00440   map<string,PetscViewer> viewers;
00441 
00442 private:
00443   PetscLogDouble start_time;    // this is used in the wall-clock-time backup code
00444 
00445   int event_step,               
00446     event_velocity,             
00447     event_energy,               
00448     event_mass,                 
00449     event_age,                  
00450     event_beddef,               
00451     event_output,               
00452     event_snapshots,            
00453     event_backups;              
00454 };
00455 
00456 #endif /* __iceModel_hh */
00457