tBedrockColumn.cc - 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
---
tBedrockColumn.cc (2869B)
---
1 /* Copyright (C) 2019 PISM Authors
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
20 #include <cassert>
21
22 #include "BedrockColumn.hh"
23
24 #include "pism/util/ConfigInterface.hh"
25
26 namespace pism {
27 namespace energy {
28
29 BedrockColumn::BedrockColumn(const std::string& prefix,
30 const Config& config, double dz, unsigned int M)
31 : m_dz(dz), m_M(M), m_system(M, prefix) {
32
33 assert(M > 1);
34
35 const double
36 rho = config.get_number("energy.bedrock_thermal.density"),
37 c = config.get_number("energy.bedrock_thermal.specific_heat_capacity");
38
39 m_k = config.get_number("energy.bedrock_thermal.conductivity");
40 m_D = m_k / (rho * c);
41 }
42
43 BedrockColumn::~BedrockColumn() {
44 // empty
45 }
46
47 /*!
48 * Advance the heat equation in time.
49 *
50 * @param[in] dt time step length
51 * @param[in] Q_bottom heat flux into the column through the bottom surface
52 * @param[in] T_top temperature at the top surface
53 * @param[in] T_old current temperature in the column
54 * @param[out] T_new output
55 *
56 * Note: T_old and T_new may point to the same location.
57 */
58 void BedrockColumn::solve(double dt, double Q_bottom, double T_top,
59 const double *T_old, double *T_new) {
60
61 double R = m_D * dt / (m_dz * m_dz);
62 double G = -Q_bottom / m_k;
63
64 m_system.L(0) = 0.0; // not used
65 m_system.D(0) = 1.0 + 2.0 * R;
66 m_system.U(0) = -2.0 * R;
67 m_system.RHS(0) = T_old[0] - 2.0 * G * m_dz * R;
68
69 unsigned int N = m_M - 1;
70
71 for (unsigned int k = 1; k < N; ++k) {
72 m_system.L(k) = -R;
73 m_system.D(k) = 1.0 + 2.0 * R;
74 m_system.U(k) = -R;
75 m_system.RHS(k) = T_old[k];
76 }
77
78 m_system.L(N) = 0.0;
79 m_system.D(N) = 1.0;
80 m_system.U(N) = 0.0; // not used
81 m_system.RHS(N) = T_top;
82
83 m_system.solve(m_M, T_new);
84 }
85
86 /*!
87 * This version of `solve()` is easier to use in Python.
88 */
89 void BedrockColumn::solve(double dt, double Q_bottom, double T_top,
90 const std::vector<double> &T_old,
91 std::vector<double> &result) {
92 result.resize(m_M);
93 solve(dt, Q_bottom, T_top, T_old.data(), result.data());
94 }
95
96
97 } // end of namespace energy
98 } // end of namespace pism