EnthalpyConverter.hh (4706B)
1 // Copyright (C) 2009-2011, 2013, 2014, 2015, 2016 Andreas Aschwanden, Ed Bueler and Constantine Khroulev 2 // 3 // This file is part of PISM. 4 // 5 // PISM is free software; you can redistribute it and/or modify it under the 6 // terms of the GNU General Public License as published by the Free Software 7 // Foundation; either version 3 of the License, or (at your option) any later 8 // version. 9 // 10 // PISM is distributed in the hope that it will be useful, but WITHOUT ANY 11 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS 12 // FOR A PARTICULAR PURPOSE. See the GNU General Public License for more 13 // details. 14 // 15 // You should have received a copy of the GNU General Public License 16 // along with PISM; if not, write to the Free Software 17 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA 18 19 #ifndef __enthalpyConverter_hh 20 #define __enthalpyConverter_hh 21 22 #include <vector> 23 #include <memory> 24 25 namespace pism { 26 27 class Config; 28 29 //! Converts between specific enthalpy and temperature or liquid content. 30 /*! 31 Use this way, for example within IceModel with Config config member: 32 \code 33 #include "enthalpyConverter.hh" 34 35 EnthalpyConverter EC(&config); // runs constructor; do after initialization of Config config 36 ... 37 for (...) { 38 ... 39 E_s = EC.enthalpy_cts(p); 40 ... etc ... 41 } 42 \endcode 43 44 The three methods that get the enthalpy from temperatures and liquid 45 fractions, namely enthalpy(), enthalpy_permissive() are more strict 46 about error checking. They throw RuntimeError if their arguments are 47 invalid. 48 49 This class is documented by [\ref AschwandenBuelerKhroulevBlatter]. 50 */ 51 class EnthalpyConverter { 52 public: 53 EnthalpyConverter(const Config &config); 54 virtual ~EnthalpyConverter(); 55 56 typedef std::shared_ptr<EnthalpyConverter> Ptr; 57 58 bool is_temperate(double E, double P) const; 59 bool is_temperate_relaxed(double E, double P) const; 60 61 double temperature(double E, double P) const; 62 double melting_temperature(double P) const; 63 double pressure_adjusted_temperature(double E, double P) const; 64 65 double water_fraction(double E, double P) const; 66 67 double enthalpy(double T, double omega, double P) const; 68 double enthalpy_cts(double P) const; 69 double enthalpy_liquid(double P) const; 70 double enthalpy_permissive(double T, double omega, double P) const; 71 72 double c() const; 73 double L(double T_m) const; 74 75 double pressure(double depth) const; 76 void pressure(const std::vector<double> &depth, 77 unsigned int ks, std::vector<double> &result) const; 78 protected: 79 void validate_E_P(double E, double P) const; 80 void validate_T_omega_P(double T, double omega, double P) const; 81 82 double temperature_cold(double E) const; 83 double enthalpy_cold(double T) const; 84 85 //! melting temperature of pure water at atmospheric pressure 86 double m_T_melting; 87 //! latent heat of fusion of water at atmospheric pressure 88 double m_L; 89 //! specific heat capacity of ice 90 double m_c_i; 91 //! specific heat capacity of pure water 92 double m_c_w; 93 //! density of ice 94 double m_rho_i; 95 //! acceleration due to gravity 96 double m_g; 97 //! atmospheric pressure 98 double m_p_air; 99 //! beta in the Clausius-Clapeyron relation (@f$ \diff{T_m}{p} = - \beta @f$). 100 double m_beta; 101 //! temperature tolerance used in `is_temperate()` in cold ice mode 102 double m_T_tolerance; 103 //! reference temperature in the definition of ice enthalpy 104 double m_T_0; 105 //! @brief if cold ice methods are selected, use `is_temperate()` 106 //! check based on temperature, not enthalpy 107 bool m_do_cold_ice_methods; 108 }; 109 110 111 //! An EnthalpyConverter for use in verification tests. 112 113 /*! Treats ice at any temperature below 10^6 Kelvin as cold (= zero liquid fraction). 114 115 The pressure dependence of the pressure-melting temperature is neglected.c; 116 117 Note: Any instance of FlowLaw uses an EnthalpyConverter; this is 118 the one used in cold mode verification code. 119 120 121 This is the special enthalpy converter that is used in 122 temperature-based verification tests only. 123 124 In these tests ice temperatures in an exact solution may exceed the 125 pressure-melting temperature, but we still want to pretend that this 126 ice is "cold" to ensure that the map from enthalpy to temperature is 127 one-to-one. (Normally enthalpy is mapped to the (temperature, water 128 fraction) pair; here water fraction is zero, so enthalpy <--> 129 (temperature, 0.) 130 131 So, I had to pick a threshold (melting) temperature that is above 132 all ice temperatures. 10^6K was chosen. 133 */ 134 class ColdEnthalpyConverter : public EnthalpyConverter { 135 public: 136 ColdEnthalpyConverter(const Config &config); 137 virtual ~ColdEnthalpyConverter(); 138 }; 139 140 } // end of namespace pism 141 142 #endif // __enthalpyConverter_hh 143