-
Notifications
You must be signed in to change notification settings - Fork 6
Expand file tree
/
Copy pathGahmSolver.cpp
More file actions
190 lines (173 loc) · 5.65 KB
/
Copy pathGahmSolver.cpp
File metadata and controls
190 lines (173 loc) · 5.65 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
// GNU General Public License v3.0
//
// This file is part of the GAHM model (https://github.com/adcirc/gahm).
// Copyright (c) 2023 ADCIRC Development Group.
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, version 3.
//
// This program is distributed in the hope that it will be useful, but
// WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
// General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.
//
// Author: Zach Cobell
// Contact: zcobell@thewaterinstitute.org
//
#include "GahmSolver.h"
#include <algorithm>
#include <cassert>
#include <cmath>
#include <cstddef>
#include <limits>
#include <string>
#include "boost/math/policies/error_handling.hpp"
#include "gahm/GahmEquations.h"
#include "physical/Atmospheric.h"
#include "physical/Earth.h"
namespace Gahm::Solver {
/**
* Construct a GAHM solver object
* @param isotach_radius Radius of the current isotach
* @param isotach_speed Speed of the current isotach
* @param vmax Maximum wind speed
* @param p_center Minimum pressure for the storm
* @param p_background Background atmospheric pressure
* @param latitude latitude of storm center
*/
GahmSolver::GahmSolver(double isotach_radius, double isotach_speed, double vmax,
double p_center, double p_background, double latitude)
: m_isotachRadius(isotach_radius),
m_isotachSpeed(isotach_speed),
m_vmax(vmax),
m_pc(p_center),
m_pbk(p_background),
m_latitude(latitude),
m_fc(Gahm::Physical::Earth::coriolis(m_latitude)),
m_rmax_guess(
GahmSolver::estimateRmax(m_pbk - m_pc, latitude, isotach_radius)),
m_rmax(std::numeric_limits<double>::max()),
m_bg(Gahm::Physical::Atmospheric::calcHollandB(m_vmax, m_pc, m_pbk)),
m_bg_tol(1e-9),
m_phi(1.0),
m_max_it(200),
m_it(0),
m_solver(m_isotachRadius, m_isotachSpeed, m_vmax, m_fc, m_bg) {}
/**
* Runs the solver and stores the solution internally
*/
void GahmSolver::solve() {
const auto guess = m_rmax_guess;
for (size_t i = 0; i < m_max_it; ++i) {
auto new_rmax = m_solver.solve(1.0, m_isotachRadius, guess);
if (!std::isnan(new_rmax) && !std::isinf(new_rmax) &&
new_rmax != std::numeric_limits<double>::max()) {
m_rmax = new_rmax;
}
assert(new_rmax > 0.0);
m_phi = GahmEquations::phi(m_vmax, m_rmax, m_bg, m_fc);
m_bg = GahmEquations::gahm_b(m_vmax, m_rmax, m_pc, m_pbk, m_fc, m_phi);
if (std::abs(m_bg - m_solver.gahm_b()) < m_bg_tol) {
m_it = i;
break;
}
if (std::isnan(new_rmax) || std::isinf(new_rmax) || std::isnan(m_bg) ||
std::isinf(m_bg) || std::isnan(m_phi) || std::isinf(m_phi)) {
throw boost::math::evaluation_error(
std::string("Solution did not converge."));
}
m_solver.setGahmB(m_bg);
}
assert(this->rmax() > 0.0);
assert(this->gahm_b() > 0.0);
}
/**
* Returns the isotach radius
* @return isotach radius
*/
auto GahmSolver::isotachRadius() const -> double { return m_isotachRadius; }
/**
* Returns the isotach speed
* @return isotach speed
*/
auto GahmSolver::isotachSpeed() const -> double { return m_isotachSpeed; }
/**
* Returns the storm latitude
* @return storm latitude
*/
auto GahmSolver::latitude() const -> double { return m_latitude; }
/**
* Returns the central pressure
* @return central pressure
*/
auto GahmSolver::p_center() const -> double { return m_pc; }
/**
* Returns the background pressure
* @return background pressure
*/
auto GahmSolver::p_background() const -> double { return m_pbk; }
/**
* Returns the coriolis force
* @return coriolis force
*/
auto GahmSolver::f_coriolis() const -> double { return m_fc; }
/**
* Returns the maximum wind speed
* @return maximum wind speed
*/
auto GahmSolver::vmax() const -> double { return m_vmax; }
/**
* Returns the solution to the radius to maximum winds. Note that when the
* solver has not been run, the solution is std::numeric_limits<double>::max()
* @return radius to maximum winds
*/
auto GahmSolver::rmax() const -> double {
if (m_it == 0) {
throw boost::math::evaluation_error(
"GAHM Solver ERROR: Solver has not run");
}
return m_rmax;
}
/**
* Returns the solution to the GAHM Holland B parameter. Note that when the
* solver has not run, the solution is the standard Holland B
* @return GAHM holland B
*/
auto GahmSolver::gahm_b() const -> double {
if (m_it == 0) {
throw boost::math::evaluation_error(
"GAHM Solver ERROR: Solver has not run");
}
return m_bg;
}
/**
* Returns the number of iterations used in the GAHM Holland B solver
* @return
*/
auto GahmSolver::it() const -> size_t { return m_it; }
/*
* Returns the solution to the GAHM phi parameter. Note that when the
* solver has not run, the solution is 1.0
*/
auto GahmSolver::phi() const -> double { return m_phi; }
/**
* Estimates the rmax. Used as the initial guess for the solver
* @param dp pressure deficit in Pascals
* @param lat latitude of the storm center
* @param isorad isotach radius that we are solving for (upper bound)
* @return estimate of rmax
*/
auto GahmSolver::estimateRmax(const double dp, const double lat,
const double isorad) -> double {
assert(dp >= 0.0);
assert(isorad > 0.0);
auto r1 =
std::exp(3.015 - 6.291e-5 * std::pow(dp / 100.0, 2.0) + 0.337 * lat);
auto r2 = 0.99 * isorad;
return std::min(r1, r2);
}
} // namespace Gahm::Solver