PISM, A Parallel Ice Sheet Model  stable v2.1-1-g6902d5502 committed by Ed Bueler on 2023-12-20 08:38:27 -0800
TemperatureModel_Verification.cc
Go to the documentation of this file.
1 /* Copyright (C) 2016, 2017, 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 "pism/verification/TemperatureModel_Verification.hh"
21 
22 #include "pism/util/error_handling.hh"
23 #include "pism/verification/tests/exactTestsFG.hh"
24 #include "pism/verification/tests/exactTestK.h"
25 #include "pism/verification/tests/exactTestO.h"
26 #include "pism/energy/utilities.hh"
27 #include "pism/util/Time.hh"
28 #include "pism/util/Context.hh"
29 
30 namespace pism {
31 namespace energy {
32 
33 static const double ApforG = 200; // m
34 static const double LforFG = 750000; // m
35 static const double ST = 1.67e-5;
36 static const double Tmin = 223.15; // K
37 
39  std::shared_ptr<const Grid> grid,
40  std::shared_ptr<const stressbalance::StressBalance> stress_balance, int testname,
41  bool bedrock_is_ice)
42  : TemperatureModel(grid, stress_balance),
43  m_testname(testname),
44  m_bedrock_is_ice(bedrock_is_ice) {
45  // empty
46 }
47 
49  const array::Scalar &ice_thickness,
50  const array::Scalar &surface_temperature,
51  const array::Scalar &climatic_mass_balance,
52  const array::Scalar &basal_heat_flux) {
53 
54  // ignore provided basal melt rate
55  (void) basal_melt_rate;
56 
58 
59  switch (m_testname) {
60  case 'F':
61  case 'G':
62  initTestFG();
63  break;
64  case 'K':
65  case 'O':
66  initTestsKO();
67  break;
68  default:
69  TemperatureModel::initialize_impl(m_basal_melt_rate, ice_thickness, surface_temperature,
70  climatic_mass_balance, basal_heat_flux);
71  }
72 
75 
76  // this will update ghosts of m_ice_enthalpy
78 }
79 
81 
83 
84  const double time = m_testname == 'F' ? 0.0 : this->time().current();
85  const double A = m_testname == 'F' ? 0.0 : ApforG;
86 
87  for (auto p = m_grid->points(); p; p.next()) {
88  const int i = p.i(), j = p.j();
89 
90  // avoid singularity at origin
91  const double r = std::max(grid::radius(*m_grid, i, j), 1.0);
92 
93  if (r > LforFG - 1.0) { // if (essentially) outside of sheet
94  m_ice_temperature.set_column(i, j, Tmin + ST * r);
95  } else {
96  TestFGParameters P = exactFG(time, r, m_grid->z(), A);
97  m_ice_temperature.set_column(i, j, P.T.data());
98  }
99  }
100 }
101 
103 
104  const unsigned int Mz = m_grid->Mz();
105 
106  std::vector<double> temperature(Mz);
107 
108  const double time = this->time().current();
109 
110  // evaluate exact solution in a column; all columns are the same
111  for (unsigned int k = 0; k < Mz; k++) {
112  if (m_testname == 'K') {
113  temperature[k] = exactK(time, m_grid->z(k), m_bedrock_is_ice ? 1 : 0).T;
114  } else {
115  temperature[k] = exactO(m_grid->z(k)).TT;
116  }
117  }
118 
119  // fill m_ice_temperature
121 
122  ParallelSection loop(m_grid->com);
123  try {
124  for (auto p = m_grid->points(); p; p.next()) {
125  m_ice_temperature.set_column(p.i(), p.j(), temperature.data());
126  }
127  } catch (...) {
128  loop.failed();
129  }
130  loop.check();
131 }
132 
133 
134 } // end of namespace energy
135 } // end of namespace pism
const Time & time() const
Definition: Component.cc:109
const std::shared_ptr< const Grid > m_grid
grid used by this component
Definition: Component.hh:156
void failed()
Indicates a failure of a parallel section.
double current() const
Current time, in seconds.
Definition: Time.cc:465
Makes sure that we call begin_access() and end_access() for all accessed array::Arrays.
Definition: Array.hh:65
void set_column(int i, int j, double c)
Set all values of scalar quantity to given a single value in a particular column.
Definition: Array3D.cc:49
void set(double c)
Result: v[j] <- c for all j.
Definition: Array.cc:707
void update_ghosts()
Updates ghost points.
Definition: Array.cc:693
const array::Scalar & basal_melt_rate() const
Basal melt rate in grounded areas. (It is set to zero elsewhere.)
Definition: EnergyModel.cc:307
array::Array3D m_ice_enthalpy
Definition: EnergyModel.hh:141
array::Scalar m_basal_melt_rate
Definition: EnergyModel.hh:143
TemperatureModel_Verification(std::shared_ptr< const Grid > grid, std::shared_ptr< const stressbalance::StressBalance > stress_balance, int testname, bool bedrock_is_ice)
void initialize_impl(const array::Scalar &basal_melt_rate, const array::Scalar &ice_thickness, const array::Scalar &surface_temperature, const array::Scalar &climatic_mass_balance, const array::Scalar &basal_heat_flux)
const array::Array3D & temperature() const
void initialize_impl(const array::Scalar &basal_melt_rate, const array::Scalar &ice_thickness, const array::Scalar &surface_temperature, const array::Scalar &climatic_mass_balance, const array::Scalar &basal_heat_flux)
struct TestKParameters exactK(double t, double z, int bedrock_is_ice)
Definition: exactTestK.c:129
struct TestOParameters exactO(double z)
Definition: exactTestO.c:112
double max(const array::Scalar &input)
Finds maximum over all the values in an array::Scalar object. Ignores ghosts.
Definition: Scalar.cc:165
static const double ApforG
static const double ST
static const double LforFG
void compute_enthalpy_cold(const array::Array3D &temperature, const array::Scalar &ice_thickness, array::Array3D &result)
Compute ice enthalpy from temperature temperature by assuming the ice has zero liquid fraction.
Definition: utilities.cc:48
static const double Tmin
double radius(const Grid &grid, int i, int j)
Returns the distance from the point (i,j) to the origin.
Definition: Grid.cc:1407
static const double k
Definition: exactTestP.cc:42
TestFGParameters exactFG(double t, double r, const std::vector< double > &z, double Cp)
Definition: exactTestsFG.cc:38
std::vector< double > T
Definition: exactTestsFG.hh:51