Loading [MathJax]/extensions/tex2jax.js
PISM, A Parallel Ice Sheet Model 2.2.2-d6b3a29ca committed by Constantine Khrulev on 2025-03-28
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
BedrockColumn.cc
Go to the documentation of this file.
1/* Copyright (C) 2019, 2023 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 "pism/energy/BedrockColumn.hh"
23#include "pism/util/ConfigInterface.hh"
24
25namespace pism {
26namespace energy {
27
28BedrockColumn::BedrockColumn(const std::string& prefix,
29 const Config& config, double dz, unsigned int M)
30 : m_dz(dz), m_M(M), m_system(M, prefix) {
31
32 assert(M > 1);
33
34 const double
35 rho = config.get_number("energy.bedrock_thermal.density"),
36 c = config.get_number("energy.bedrock_thermal.specific_heat_capacity");
37
38 m_k = config.get_number("energy.bedrock_thermal.conductivity");
39 m_D = m_k / (rho * c);
40}
41
42/*!
43 * Advance the heat equation in time.
44 *
45 * @param[in] dt time step length
46 * @param[in] Q_bottom heat flux into the column through the bottom surface
47 * @param[in] T_top temperature at the top surface
48 * @param[in] T_old current temperature in the column
49 * @param[out] T_new output
50 *
51 * Note: T_old and T_new may point to the same location.
52 */
53void BedrockColumn::solve(double dt, double Q_bottom, double T_top,
54 const double *T_old, double *T_new) {
55
56 double R = m_D * dt / (m_dz * m_dz);
57 double G = -Q_bottom / m_k;
58
59 m_system.L(0) = 0.0; // not used
60 m_system.D(0) = 1.0 + 2.0 * R;
61 m_system.U(0) = -2.0 * R;
62 m_system.RHS(0) = T_old[0] - 2.0 * G * m_dz * R;
63
64 unsigned int N = m_M - 1;
65
66 for (unsigned int k = 1; k < N; ++k) {
67 m_system.L(k) = -R;
68 m_system.D(k) = 1.0 + 2.0 * R;
69 m_system.U(k) = -R;
70 m_system.RHS(k) = T_old[k];
71 }
72
73 m_system.L(N) = 0.0;
74 m_system.D(N) = 1.0;
75 m_system.U(N) = 0.0; // not used
76 m_system.RHS(N) = T_top;
77
78 m_system.solve(m_M, T_new);
79}
80
81/*!
82 * This version of `solve()` is easier to use in Python.
83 */
84void BedrockColumn::solve(double dt, double Q_bottom, double T_top,
85 const std::vector<double> &T_old,
86 std::vector<double> &result) {
87 result.resize(m_M);
88 solve(dt, Q_bottom, T_top, T_old.data(), result.data());
89}
90
91
92} // end of namespace energy
93} // end of namespace pism
double get_number(const std::string &name, UseFlag flag=REMEMBER_THIS_USE) const
A class for storing and accessing PISM configuration flags and parameters.
double & D(size_t i)
double & L(size_t i)
double & RHS(size_t i)
void solve(unsigned int system_size, std::vector< double > &result)
double & U(size_t i)
BedrockColumn(const std::string &prefix, const Config &config, double dz, unsigned int M)
void solve(double dt, double Q_bottom, double T_top, const double *T_old, double *result)
const double G
Definition exactTestK.c:40
#define rho
Definition exactTestM.c:35
static const double k
Definition exactTestP.cc:42