ARTS 2.5.4 (git: 31ce4f0e)
m_xsec_fit.cc
Go to the documentation of this file.
1/* Copyright (C) 2017 Oliver Lemke <oliver.lemke@uni-hamburg.de> and
2 Stefan Buehler <stefan.buehler@uni-hamburg.de>.
3
4 This program is free software; you can redistribute it and/or modify it
5 under the terms of the GNU General Public License as published by the
6 Free Software Foundation; either version 2, or (at your option) any
7 later version.
8
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
13
14 You should have received a copy of the GNU General Public License
15 along with this program; if not, write to the Free Software
16 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
17 USA. */
18
28#include "arts.h"
29#include "global_data.h"
30#include "xsec_fit.h"
31#include "jacobian.h"
32#include "m_xml.h"
33#include "messages.h"
34#include "physics_funcs.h"
35
36/* Workspace method: Doxygen documentation will be auto-generated */
37void ReadXsecData(ArrayOfXsecRecord& xsec_fit_data,
38 const ArrayOfArrayOfSpeciesTag& abs_species,
39 const String& basename,
40 const Verbosity& verbosity) {
41 // Build a set of species indices. Duplicates are ignored.
42 std::set<Species::Species> unique_species;
43 for (auto& asp : abs_species) {
44 for (auto& sp : asp) {
45 if (sp.Type() == Species::TagType::XsecFit) {
46 unique_species.insert(sp.Spec());
47 }
48 }
49 }
50
51 String tmpbasename = basename;
52 if (basename.length() && basename[basename.length() - 1] != '/') {
53 tmpbasename += '.';
54 }
55
56 // Read xsec data for all active species and collect them in xsec_fit_data
57 xsec_fit_data.clear();
58 for (auto& species_name : unique_species) {
59 XsecRecord xsec_coeffs;
60 const String filename{tmpbasename +
61 String(Species::toShortName(species_name)) + ".xml"};
62
63 try {
64 ReadXML(xsec_coeffs, "", filename, "", verbosity);
65
66 xsec_fit_data.push_back(xsec_coeffs);
67 } catch (const std::exception& e) {
69 "Error reading coefficients file:\n", filename, "\n", e.what());
70 }
71 }
72}
73
74/* Workspace method: Doxygen documentation will be auto-generated */
75void propmat_clearskyAddXsecFit( // WS Output:
76 PropagationMatrix& propmat_clearsky,
77 ArrayOfPropagationMatrix& dpropmat_clearsky_dx,
78 // WS Input:
79 const ArrayOfArrayOfSpeciesTag& abs_species,
80 const ArrayOfSpeciesTag& select_abs_species,
81 const ArrayOfRetrievalQuantity& jacobian_quantities,
82 const Vector& f_grid,
83 const Numeric& rtp_pressure,
84 const Numeric& rtp_temperature,
85 const Vector& rtp_vmr,
86 const ArrayOfXsecRecord& xsec_fit_data,
87 const Numeric& force_p,
88 const Numeric& force_t,
89 // Verbosity object:
90 const Verbosity& verbosity) {
92
93 // Forward simulations and their error handling
94 ARTS_USER_ERROR_IF(rtp_vmr.nelem() not_eq abs_species.nelem(),
95 "Mismatch dimensions on species and VMR inputs");
97 propmat_clearsky.NumberOfFrequencies() not_eq f_grid.nelem(),
98 "Mismatch dimensions on internal matrices of xsec and frequency");
99
100 // Derivatives and their error handling
101 if (dpropmat_clearsky_dx.nelem()) {
103 dpropmat_clearsky_dx.nelem() not_eq jacobian_quantities.nelem(),
104 "Mismatch dimensions on xsec derivatives and Jacobian grids");
106 std::any_of(dpropmat_clearsky_dx.cbegin(),
107 dpropmat_clearsky_dx.cend(),
108 [&f_grid](auto& x) {
109 return x.NumberOfFrequencies() not_eq f_grid.nelem();
110 }),
111 "Mismatch dimensions on internal matrices of xsec derivatives and frequency");
112 }
113
114 // Jacobian overhead START
115 /* NOTE: The calculations below are inefficient and could
116 be made much better by using interp in Extract to
117 return the derivatives as well. */
118 // Jacobian vectors START
119 Vector dxsec_temp_dT;
120 Vector dxsec_temp_dF;
121 Vector dfreq;
122 // Jacobian vectors END
123 const bool do_jac = supports_hitran_xsec(jacobian_quantities);
124 const bool do_freq_jac = do_frequency_jacobian(jacobian_quantities);
125 const bool do_temp_jac = do_temperature_jacobian(jacobian_quantities);
126 const Numeric df = frequency_perturbation(jacobian_quantities);
127 const Numeric dt = temperature_perturbation(jacobian_quantities);
128 if (do_freq_jac) {
129 dfreq.resize(f_grid.nelem());
130 dfreq = f_grid;
131 dfreq += df;
132 dxsec_temp_dF.resize(f_grid.nelem());
133 }
134 if (do_temp_jac) {
135 dxsec_temp_dT.resize(f_grid.nelem());
136 }
137 // Jacobian overhead END
138
139 // Useful if there is no Jacobian to calculate
140 ArrayOfMatrix empty;
141
142 // Allocate a vector with dimension frequencies for constructing our
143 // cross-sections before adding them (more efficient to allocate this here
144 // outside of the loops)
145 Vector xsec_temp(f_grid.nelem(), 0.);
146
147 ArrayOfString fail_msg;
148 bool do_abort = false;
149 // Loop over Xsec data sets.
150 // Index ii loops through the outer array (different tag groups),
151 // index s through the inner array (different tags within each goup).
152 for (Index i = 0; i < abs_species.nelem(); i++) {
153 if (select_abs_species.nelem() and abs_species[i] not_eq select_abs_species)
154 continue;
155
156 for (Index s = 0; s < abs_species[i].nelem(); s++) {
157 const SpeciesTag& this_species = abs_species[i][s];
158
159 // Check if this is a HITRAN cross section tag
160 if (this_species.Type() != Species::TagType::XsecFit) continue;
161
162 Index this_xdata_index =
163 hitran_xsec_get_index(xsec_fit_data, this_species.Spec());
164 ARTS_USER_ERROR_IF(this_xdata_index < 0,
165 "Cross-section species ",
166 this_species.Name(),
167 " not found in *xsec_fit_data*.")
168 const XsecRecord& this_xdata = xsec_fit_data[this_xdata_index];
169 // ArrayOfMatrix& this_dxsec = do_jac ? dpropmat_clearsky_dx[i] : empty;
170
171 if (do_abort) continue;
172 const Numeric current_p = force_p < 0 ? rtp_pressure : force_p;
173 const Numeric current_t = force_t < 0 ? rtp_temperature : force_t;
174
175 // Get the absorption cross sections from the HITRAN data:
176 this_xdata.Extract(xsec_temp, f_grid, current_p, current_t, verbosity);
177 if (do_freq_jac) {
178 this_xdata.Extract(
179 dxsec_temp_dF, dfreq, current_p, current_t, verbosity);
180 }
181 if (do_temp_jac) {
182 this_xdata.Extract(
183 dxsec_temp_dT, f_grid, current_p, current_t + dt, verbosity);
184 }
185 }
186
187 // Add to result variable:
188 Numeric nd = number_density(rtp_pressure, rtp_temperature);
189 if (!do_jac) {
190 xsec_temp *= nd * rtp_vmr[i];
191 propmat_clearsky.Kjj() += xsec_temp;
192 } else {
193 Numeric dnd_dt = dnumber_density_dt(rtp_pressure, rtp_temperature);
194 for (Index f = 0; f < f_grid.nelem(); f++) {
195 propmat_clearsky.Kjj()[f] += xsec_temp[f] * nd * rtp_vmr[i];
196 for (Index iq = 0; iq < jacobian_quantities.nelem(); iq++) {
197 const auto& deriv = jacobian_quantities[iq];
198
199 if (!deriv.propmattype()) continue;
200
201 if (is_frequency_parameter(deriv)) {
202 dpropmat_clearsky_dx[iq].Kjj()[f] +=
203 (dxsec_temp_dF[f] - xsec_temp[f]) * nd * rtp_vmr[i] / df;
204 } else if (deriv == Jacobian::Special::ArrayOfSpeciesTagVMR ||
205 deriv == Jacobian::Line::VMR) {
206 if (species_match(deriv, abs_species[i])) {
207 dpropmat_clearsky_dx[iq].Kjj()[f] +=
208 xsec_temp[f] * nd * rtp_vmr[i];
209 }
210 } else if (deriv == Jacobian::Atm::Temperature) {
211 dpropmat_clearsky_dx[iq].Kjj()[f] +=
212 ((dxsec_temp_dT[f] - xsec_temp[f]) / dt * nd +
213 xsec_temp[f] * dnd_dt) *
214 rtp_vmr[i];
215 }
216 }
217 }
218 }
219 }
220}
221
222/* Workspace method: Doxygen documentation will be auto-generated */
224 ArrayOfMatrix& abs_xsec_per_species,
225 ArrayOfArrayOfMatrix& /* dabs_xsec_per_species_dx */,
226 // WS Input:
227 const ArrayOfArrayOfSpeciesTag& abs_species,
228 const ArrayOfRetrievalQuantity& /* jacobian_quantities */,
229 const ArrayOfIndex& abs_species_active,
230 const Vector& f_grid,
231 const Vector& abs_p,
232 const Vector& abs_t,
233 const ArrayOfXsecRecord& xsec_fit_data,
234 const Numeric& force_p,
235 const Numeric& force_t,
236 // Verbosity object:
237 const Verbosity& verbosity) {
239
240 {
241 // Check that all parameters that should have the number of tag
242 // groups as a dimension are consistent:
243 const Index n_tgs = abs_species.nelem();
244 const Index n_xsec = abs_xsec_per_species.nelem();
245
247 n_tgs != n_xsec,
248 "The following variables must all have the same dimension:\n"
249 "abs_species: ",
250 n_tgs,
251 "\n"
252 "abs_xsec_per_species: ",
253 n_xsec)
254 }
255
256 // Useful if there is no Jacobian to calculate
257 ArrayOfMatrix empty;
258
259 {
260 // Check that all parameters that should have the the dimension of p_grid
261 // are consistent:
262 const Index n_p = abs_p.nelem();
263 const Index n_t = abs_t.nelem();
264
266 n_p != n_t,
267 "The following variables must all have the same dimension:\n"
268 "abs_p: ",
269 n_p,
270 "\n"
271 "abs_t: ",
272 n_t)
273 }
274
275 // Allocate a vector with dimension frequencies for constructing our
276 // cross-sections before adding them (more efficient to allocate this here
277 // outside of the loops)
278 Vector xsec_temp(f_grid.nelem(), 0.);
279
280 ArrayOfString fail_msg;
281 bool do_abort = false;
282
283 // Loop over Xsec data sets.
284 // Index ii loops through the outer array (different tag groups),
285 // index s through the inner array (different tags within each goup).
286 for (Index ii = 0; ii < abs_species_active.nelem(); ii++) {
287 const Index i = abs_species_active[ii];
288
289 for (Index s = 0; s < abs_species[i].nelem(); s++) {
290 const SpeciesTag& this_species = abs_species[i][s];
291
292 // Check if this is a HITRAN cross section tag
293 if (this_species.Type() != Species::TagType::XsecFit) continue;
294
295 Index this_xdata_index =
296 hitran_xsec_get_index(xsec_fit_data, this_species.Spec());
297 ARTS_USER_ERROR_IF(this_xdata_index < 0,
298 "Cross-section species ",
299 this_species.Name(),
300 " not found in *xsec_fit_data*.")
301 const XsecRecord& this_xdata = xsec_fit_data[this_xdata_index];
302 Matrix& this_xsec = abs_xsec_per_species[i];
303
304 // Loop over pressure:
305#pragma omp parallel for if (!arts_omp_in_parallel() && abs_p.nelem() >= 1) \
306 firstprivate(xsec_temp)
307 for (Index ip = 0; ip < abs_p.nelem(); ip++) {
308 if (do_abort) continue;
309 const Numeric current_p = force_p < 0 ? abs_p[ip] : force_p;
310 const Numeric current_t = force_t < 0 ? abs_t[ip] : force_t;
311
312 // Get the absorption cross sections from the HITRAN data:
313 try {
314 this_xdata.Extract(xsec_temp,
315 f_grid,
316 current_p,
317 current_t,
318 verbosity);
319 } catch (runtime_error& e) {
320 ostringstream os;
321 os << "Problem with HITRAN cross section species "
322 << this_species.Name() << " at pressure level " << ip << " ("
323 << abs_p[ip] / 100. << " hPa):\n"
324 << e.what() << "\n";
325#pragma omp critical(abs_xsec_per_speciesAddXsecFit)
326 {
327 do_abort = true;
328 fail_msg.push_back(os.str());
329 }
330 }
331 this_xsec(joker, ip) += xsec_temp;
332 }
333 }
334 }
335
336 if (do_abort) {
337 ARTS_USER_ERROR(fail_msg);
338 }
339}
The global header file for ARTS.
This can be used to make arrays out of anything.
Definition: array.h:48
Index nelem() const ARTS_NOEXCEPT
Definition: array.h:92
Index nelem() const noexcept
Returns the number of elements.
Definition: matpackI.h:536
The Matrix class.
Definition: matpackI.h:1261
Index NumberOfFrequencies() const
The number of frequencies of the propagation matrix.
VectorView Kjj(const Index iz=0, const Index ia=0)
Vector view to diagonal elements.
The Vector class.
Definition: matpackI.h:899
void resize(Index n)
Resize function.
Definition: matpackI.cc:390
Hitran crosssection class.
Definition: xsec_fit.h:43
void Extract(VectorView result, const Vector &f_grid, Numeric pressure, Numeric temperature, const Verbosity &verbosity) const
Calculate hitran cross section data.
Definition: xsec_fit.cc:52
#define ARTS_USER_ERROR(...)
Definition: debug.h:150
#define ARTS_USER_ERROR_IF(condition,...)
Definition: debug.h:134
bool do_frequency_jacobian(const ArrayOfRetrievalQuantity &js) noexcept
Returns if the array wants a frequency derivative.
Definition: jacobian.cc:1175
bool species_match(const RetrievalQuantity &rq, const ArrayOfSpeciesTag &ast)
Returns if the Retrieval quantity is VMR derivative for all the species in the species tags.
Definition: jacobian.cc:1119
bool do_temperature_jacobian(const ArrayOfRetrievalQuantity &js) noexcept
Returns if the array wants the temperature derivative.
Definition: jacobian.cc:1150
bool supports_hitran_xsec(const ArrayOfRetrievalQuantity &js)
Returns if the array supports HITRAN cross-section derivatives.
Definition: jacobian.cc:1093
bool is_frequency_parameter(const RetrievalQuantity &t) noexcept
Returns if the Retrieval quantity is a frequency parameter in propagation matrix calculations.
Definition: jacobian.cc:1006
Numeric temperature_perturbation(const ArrayOfRetrievalQuantity &js) noexcept
Returns the temperature perturbation if it exists.
Definition: jacobian.cc:1183
Numeric frequency_perturbation(const ArrayOfRetrievalQuantity &js) noexcept
Returns the frequency perturbation if it exists.
Definition: jacobian.cc:1191
Routines for setting up the jacobian.
Workspace methods and template functions for supergeneric XML IO.
void ReadXML(T &v, const String &v_name, const String &f, const String &f_name, const Verbosity &verbosity)
WORKSPACE METHOD: ReadXML.
Definition: m_xml.h:41
void abs_xsec_per_speciesAddXsecFit(ArrayOfMatrix &abs_xsec_per_species, ArrayOfArrayOfMatrix &, const ArrayOfArrayOfSpeciesTag &abs_species, const ArrayOfRetrievalQuantity &, const ArrayOfIndex &abs_species_active, const Vector &f_grid, const Vector &abs_p, const Vector &abs_t, const ArrayOfXsecRecord &xsec_fit_data, const Numeric &force_p, const Numeric &force_t, const Verbosity &verbosity)
Definition: m_xsec_fit.cc:223
void propmat_clearskyAddXsecFit(PropagationMatrix &propmat_clearsky, ArrayOfPropagationMatrix &dpropmat_clearsky_dx, const ArrayOfArrayOfSpeciesTag &abs_species, const ArrayOfSpeciesTag &select_abs_species, const ArrayOfRetrievalQuantity &jacobian_quantities, const Vector &f_grid, const Numeric &rtp_pressure, const Numeric &rtp_temperature, const Vector &rtp_vmr, const ArrayOfXsecRecord &xsec_fit_data, const Numeric &force_p, const Numeric &force_t, const Verbosity &verbosity)
WORKSPACE METHOD: propmat_clearskyAddXsecFit.
Definition: m_xsec_fit.cc:75
void ReadXsecData(ArrayOfXsecRecord &xsec_fit_data, const ArrayOfArrayOfSpeciesTag &abs_species, const String &basename, const Verbosity &verbosity)
WORKSPACE METHOD: ReadXsecData.
Definition: m_xsec_fit.cc:37
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
const Joker joker
Declarations having to do with the four output streams.
#define CREATE_OUTS
Definition: messages.h:208
my_basic_string< char > String
The String type for ARTS.
Definition: mystring.h:216
constexpr Numeric e
Elementary charge convenience name [C].
Particulates ENUMCLASS(Line, char, VMR, Strength, Center, ShapeG0X0, ShapeG0X1, ShapeG0X2, ShapeG0X3, ShapeD0X0, ShapeD0X1, ShapeD0X2, ShapeD0X3, ShapeG2X0, ShapeG2X1, ShapeG2X2, ShapeG2X3, ShapeD2X0, ShapeD2X1, ShapeD2X2, ShapeD2X3, ShapeFVCX0, ShapeFVCX1, ShapeFVCX2, ShapeFVCX3, ShapeETAX0, ShapeETAX1, ShapeETAX2, ShapeETAX3, ShapeYX0, ShapeYX1, ShapeYX2, ShapeYX3, ShapeGX0, ShapeGX1, ShapeGX2, ShapeGX3, ShapeDVX0, ShapeDVX1, ShapeDVX2, ShapeDVX3, ECS_SCALINGX0, ECS_SCALINGX1, ECS_SCALINGX2, ECS_SCALINGX3, ECS_BETAX0, ECS_BETAX1, ECS_BETAX2, ECS_BETAX3, ECS_LAMBDAX0, ECS_LAMBDAX1, ECS_LAMBDAX2, ECS_LAMBDAX3, ECS_DCX0, ECS_DCX1, ECS_DCX2, ECS_DCX3, NLTE) static_assert(Index(Line ArrayOfSpeciesTagVMR
Definition: jacobian.h:102
Temperature
Definition: jacobian.h:59
This file contains declerations of functions of physical character.
constexpr Numeric dnumber_density_dt(Numeric p, Numeric t) noexcept
dnumber_density_dT
Definition: physics_funcs.h:84
constexpr Numeric number_density(Numeric p, Numeric t) noexcept
number_density
Definition: physics_funcs.h:70
Species::Tag SpeciesTag
Definition: species_tags.h:101
Index hitran_xsec_get_index(const ArrayOfXsecRecord &xsec_data, const Species::Species species)
Get the index in xsec_fit_data for the given species.
Definition: xsec_fit.cc:226
Methods and classes for HITRAN absorption cross section data.