tenthSystem.hh - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
HTML git clone git://src.adamsgaard.dk/pism
DIR Log
DIR Files
DIR Refs
DIR LICENSE
---
tenthSystem.hh (3758B)
---
1 // Copyright (C) 2009-2011, 2013, 2014, 2015, 2016, 2017, 2018 Andreas Aschwanden and Ed Bueler
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 __enthSystem_hh
20 #define __enthSystem_hh
21
22 #include <vector>
23
24 #include "pism/util/ColumnSystem.hh"
25 #include "pism/util/EnthalpyConverter.hh"
26
27 namespace pism {
28
29 class Config;
30 class IceModelVec3;
31
32 namespace energy {
33
34 //! Tridiagonal linear system for conservation of energy in vertical column of ice enthalpy.
35 /*!
36 See the page documenting \ref bombproofenth. The top of
37 the ice has a Dirichlet condition.
38 */
39 class enthSystemCtx : public columnSystemCtx {
40
41 public:
42 enthSystemCtx(const std::vector<double>& storage_grid,
43 const std::string &prefix,
44 double dx, double dy, double dt,
45 const Config &config,
46 const IceModelVec3 &Enth3,
47 const IceModelVec3 &u3,
48 const IceModelVec3 &v3,
49 const IceModelVec3 &w3,
50 const IceModelVec3 &strain_heating3,
51 EnthalpyConverter::Ptr EC);
52 ~enthSystemCtx();
53
54 void init(int i, int j, bool ismarginal, double ice_thickness);
55
56 double k_from_T(double T) const;
57
58 void set_surface_heat_flux(double hf);
59 void set_surface_neumann_bc(double dE);
60 void set_surface_dirichlet_bc(double E_surface);
61
62 void set_basal_dirichlet_bc(double E_basal);
63 void set_basal_heat_flux(double hf);
64 void set_basal_neumann_bc(double dE);
65
66 virtual void save_system(std::ostream &output, unsigned int M) const;
67
68 void solve(std::vector<double> &result);
69
70 double lambda() const {
71 return m_lambda;
72 }
73
74 double Enth(size_t i) const {
75 return m_Enth[i];
76 }
77
78 double Enth_s(size_t i) const {
79 return m_Enth_s[i];
80 }
81 protected:
82 // enthalpy in ice at previous time step
83 std::vector<double> m_Enth;
84 // enthalpy level for CTS; function only of pressure
85 std::vector<double> m_Enth_s;
86
87 // temporary storage for ice enthalpy at (i,j), as well as north,
88 // east, south, and west from (i,j)
89 std::vector<double> m_E_ij, m_E_n, m_E_e, m_E_s, m_E_w;
90
91 //! strain heating in the ice column
92 std::vector<double> m_strain_heating;
93
94 //! values of @f$ k \Delta t / (\rho c \Delta x^2) @f$
95 std::vector<double> m_R;
96
97 double m_ice_density, m_ice_c, m_ice_k, m_p_air,
98 m_nu, m_R_cold, m_R_temp, m_R_factor;
99
100 double m_ice_thickness,
101 m_lambda; //!< implicit FD method parameter
102 double m_D0, m_U0, m_B0; // coefficients of the first (basal) equation
103 double m_L_ks, m_D_ks, m_U_ks, m_B_ks; // coefficients of the last (surface) equation
104 bool m_marginal, m_k_depends_on_T;
105
106 bool m_exclude_horizontal_advection;
107 bool m_exclude_vertical_advection;
108 bool m_exclude_strain_heat;
109
110 const IceModelVec3 &m_Enth3, &m_strain_heating3;
111 EnthalpyConverter::Ptr m_EC; // conductivity has known dependence on T, not enthalpy
112
113 void compute_enthalpy_CTS();
114 double compute_lambda();
115
116 void assemble_R();
117 void checkReadyToSolve();
118 };
119
120 } // end of namespace energy
121 } // end of namespace pism
122
123 #endif // ifndef __enthSystem_hh