ARTS 2.5.10 (git: 2f1c442c)
m_nlte.cc
Go to the documentation of this file.
1/* Copyright (C) 2018
2 Richard Larsson
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
27#include "absorption.h"
28#include "arts.h"
29#include "auto_md.h"
30#include "lin_alg.h"
31#include "nlte.h"
32#include "xml_io.h"
33
34/* Workspace method: Doxygen documentation will be auto-generated */
37 const ArrayOfArrayOfAbsorptionLines& abs_lines_per_species,
38 const Index& global,
39 const Verbosity&) {
40 // Defined only as output not input so resizing
41 qid.resize(0);
42
43 QuantumIdentifier lower, upper;
44
45 // For all lines' upper and lower energy levels
46 for (const auto& lines: abs_lines_per_species) {
47 for (const auto& band: lines) {
48 for (Index k=0; k<band.NumLines() and (global ? (k==0) : false); k++) {
49 if (global) {
50 lower = band.quantumidentity.LowerLevel();
51 upper = band.quantumidentity.UpperLevel();
52 } else {
53 auto x = band.QuantumIdentityOfLine(k);
54 lower = x.LowerLevel();
55 upper = x.UpperLevel();
56 }
57
58 if (std::none_of(qid.begin(), qid.end(), [&](auto& x){return x == lower;}))
59 qid.push_back(lower);
60 if (std::none_of(qid.begin(), qid.end(), [&](auto& x){return x == upper;}))
61 qid.push_back(upper);
62 }
63 }
64 }
65}
66
67/* Workspace method: Doxygen documentation will be auto-generated */
69 const Numeric& scale,
70 const Verbosity&) {
71 nlte_field.value *= scale;
72}
73
74/* Workspace method: Doxygen documentation will be auto-generated */
76 Workspace& ws,
77 EnergyLevelMap& nlte_field,
78 const ArrayOfArrayOfSpeciesTag& abs_species,
79 const ArrayOfArrayOfAbsorptionLines& abs_lines_per_species,
80 const ArrayOfArrayOfGriddedField1& collision_coefficients,
81 const ArrayOfQuantumIdentifier& collision_line_identifiers,
82 const SpeciesIsotopologueRatios& isotopologue_ratios,
83 const Agenda& iy_main_agenda,
84 const Agenda& ppath_agenda,
85 const Agenda& iy_space_agenda,
86 const Agenda& iy_surface_agenda,
87 const Agenda& iy_cloudbox_agenda,
88 const Agenda& propmat_clearsky_agenda,
89 const Agenda& /*water_p_eq_agenda*/,
90 const Tensor4& vmr_field,
91 const Tensor3& t_field,
92 const Tensor3& z_field,
93 const Vector& p_grid,
94 const Index& atmosphere_dim,
95 const Vector& refellipsoid,
96 const Tensor3& surface_props_data,
97 const Index& nlte_do,
98 const Numeric& df,
99 const Numeric& convergence_limit,
100 const Index& nz,
101 const Index& nf,
102 const Index& dampened,
103 const Index& iteration_limit,
104 const Verbosity& verbosity)
105{
107
108 ARTS_USER_ERROR_IF (not nlte_do, "Must be set to do NLTE");
109 ARTS_USER_ERROR_IF (nlte_field.value.empty(),
110 "Error in NLTE field, it is empty");
111
112 Matrix line_irradiance;
113 Tensor3 line_transmission;
114
115 const Index nlevels = nlte_field.levels.nelem(), np = p_grid.nelem();
116 ARTS_USER_ERROR_IF (nlevels < 5,
117 "Must have more than a four levels");
118
119 ARTS_USER_ERROR_IF (atmosphere_dim not_eq 1,
120 "Only for 1D atmosphere");
121
122 const Index nlines = nelem(abs_lines_per_species);
123 ARTS_USER_ERROR_IF (nlevels >= nlines,
124 "Bad number of lines... overlapping lines in nlte_level_identifiers?");
125
126 // Create basic compute vectors
127 const Vector Aij = createAij(abs_lines_per_species);
128 const Vector Bij = createBij(abs_lines_per_species);
129 const Vector Bji = createBji(Bij, abs_lines_per_species);
130 Vector Cij(nlines), Cji(nlines);
131
132 ArrayOfIndex upper, lower;
134 upper, lower, abs_lines_per_species, nlte_field);
135 const Index unique = find_first_unique_in_lower(upper, lower);
136
137 // Compute arrays
138 Matrix SEE(nlevels, nlevels, 0.0);
139 Vector r(nlevels, 0.0), x(nlevels, 0.0);
140 Numeric max_change = convergence_limit + 1;
141
142 Index i = 0;
143 while (i < iteration_limit and max_change > convergence_limit) {
144 // Reset change
145 max_change = 0.0;
146
147 // Compute radiation and transmission
149 ws,
150 line_irradiance,
151 line_transmission,
152 abs_species,
153 abs_lines_per_species,
154 nlte_field,
155 vmr_field,
156 t_field,
157 z_field,
158 p_grid,
159 refellipsoid,
160 surface_props_data,
161 ppath_agenda,
162 iy_main_agenda,
163 iy_space_agenda,
164 iy_surface_agenda,
165 iy_cloudbox_agenda,
166 propmat_clearsky_agenda,
167 df,
168 nz,
169 nf,
170 1.0,
171 verbosity);
172
173 for (Index ip = 0; ip < np; ip++) {
174 r = nlte_field.value(joker, ip, 0, 0);
176 Cji,
177 abs_lines_per_species,
178 abs_species,
179 collision_coefficients,
180 collision_line_identifiers,
181 isotopologue_ratios,
182 vmr_field(joker, ip, 0, 0),
183 t_field(ip, 0, 0),
184 p_grid[ip]);
185
186
187 if (dampened)
189 SEE,
190 r,
191 Aij,
192 Bij,
193 Bji,
194 Cij,
195 Cji,
196 line_irradiance(joker, ip),
197 line_transmission(0, joker, ip),
198 upper,
199 lower);
200 else
202 Aij,
203 Bij,
204 Bji,
205 Cij,
206 Cji,
207 line_irradiance(joker, ip),
208 upper,
209 lower);
210
212 solve(nlte_field.value(joker, ip, 0, 0), SEE, x);
213
214 for (Index il = 0; il < nlevels; il++) {
215 max_change =
216 max(abs(nlte_field.value(il, ip, 0, 0) - r[il]) / r[il], max_change);
217 }
218 }
219 i++;
220 }
221
222 if (i < iteration_limit)
223 out2 << "Converged NLTE ratios (within convergence_limit) returned after "
224 << i << " iterations\n";
225 else
226 out2
227 << "No convergence of NLTE ratios (within convergence_limit) returned even after "
228 << iteration_limit << " iterations\n";
229}
230
231/* Workspace method: Doxygen documentation will be auto-generated */
233 ArrayOfArrayOfGriddedField1& collision_coefficients,
234 ArrayOfQuantumIdentifier& collision_line_identifiers,
235 const ArrayOfArrayOfSpeciesTag& abs_species,
236 const String& basename,
237 const Verbosity& verbosity) {
238 // Standard ARTS basename-assumption
239 String tmp_basename = basename, filename;
240 if (basename.length() && basename[basename.length() - 1] != '/')
241 tmp_basename += ".";
242
243 // Read the identification file
244 filename = tmp_basename + "qid.xml";
245 xml_read_from_file(filename, collision_line_identifiers, verbosity);
246 check_collision_line_identifiers(collision_line_identifiers);
247
248 // Inner array size has to be this constantly
249 const Index n = collision_line_identifiers.nelem();
250
251 // Set species dimensions and fill the array
252 collision_coefficients.resize(abs_species.nelem());
253 for (Index i = 0; i < collision_coefficients.nelem(); i++) {
255
256 // Read the file for a species and check that the size is correct of the array
257 filename = tmp_basename + String(Species::toShortName(abs_species[i].Species())) + ".xml";
258 xml_read_from_file(filename, aogf1, verbosity);
259 ARTS_USER_ERROR_IF (aogf1.nelem() not_eq n,
260 "Mismatch between collision_line_identifiers and some collision_coefficients");
261 collision_coefficients[i] = aogf1;
262 }
263}
264
265/* Workspace method: Doxygen documentation will be auto-generated */
266void nlteOff(Index& nlte_do,
267 EnergyLevelMap& nlte_field,
268 ArrayOfQuantumIdentifier& nlte_level_identifiers,
269 const Verbosity&) {
270 nlte_do = 0;
271 nlte_field = EnergyLevelMap();
272 nlte_level_identifiers.resize(0);
273}
Declarations required for the calculation of absorption coefficients.
base max(const Array< base > &x)
Max function.
Definition: array.h:145
The global header file for ARTS.
The Agenda class.
Definition: agenda_class.h:69
Index nelem() const ARTS_NOEXCEPT
Definition: array.h:92
bool empty() const noexcept
Definition: matpackIV.h:152
Numeric sum() const ARTS_NOEXCEPT
The sum of all elements of a Vector.
Definition: matpackI.cc:48
Index nelem() const noexcept
Returns the number of elements.
Definition: matpackI.h:547
The Matrix class.
Definition: matpackI.h:1285
The Tensor3 class.
Definition: matpackIII.h:352
The Tensor4 class.
Definition: matpackIV.h:435
The Vector class.
Definition: matpackI.h:910
Workspace class.
Definition: workspace_ng.h:53
void solve(VectorView w, const CovarianceMatrix &A, ConstVectorView v)
#define ARTS_USER_ERROR_IF(condition,...)
Definition: debug.h:153
Linear algebra functions.
void nlte_fieldForSingleSpeciesNonOverlappingLines(Workspace &ws, EnergyLevelMap &nlte_field, const ArrayOfArrayOfSpeciesTag &abs_species, const ArrayOfArrayOfAbsorptionLines &abs_lines_per_species, const ArrayOfArrayOfGriddedField1 &collision_coefficients, const ArrayOfQuantumIdentifier &collision_line_identifiers, const SpeciesIsotopologueRatios &isotopologue_ratios, const Agenda &iy_main_agenda, const Agenda &ppath_agenda, const Agenda &iy_space_agenda, const Agenda &iy_surface_agenda, const Agenda &iy_cloudbox_agenda, const Agenda &propmat_clearsky_agenda, const Agenda &, const Tensor4 &vmr_field, const Tensor3 &t_field, const Tensor3 &z_field, const Vector &p_grid, const Index &atmosphere_dim, const Vector &refellipsoid, const Tensor3 &surface_props_data, const Index &nlte_do, const Numeric &df, const Numeric &convergence_limit, const Index &nz, const Index &nf, const Index &dampened, const Index &iteration_limit, const Verbosity &verbosity)
WORKSPACE METHOD: nlte_fieldForSingleSpeciesNonOverlappingLines.
Definition: m_nlte.cc:75
void collision_coefficientsFromSplitFiles(ArrayOfArrayOfGriddedField1 &collision_coefficients, ArrayOfQuantumIdentifier &collision_line_identifiers, const ArrayOfArrayOfSpeciesTag &abs_species, const String &basename, const Verbosity &verbosity)
WORKSPACE METHOD: collision_coefficientsFromSplitFiles.
Definition: m_nlte.cc:232
void nlteOff(Index &nlte_do, EnergyLevelMap &nlte_field, ArrayOfQuantumIdentifier &nlte_level_identifiers, const Verbosity &)
WORKSPACE METHOD: nlteOff.
Definition: m_nlte.cc:266
void nlte_fieldRescalePopulationLevels(EnergyLevelMap &nlte_field, const Numeric &scale, const Verbosity &)
WORKSPACE METHOD: nlte_fieldRescalePopulationLevels.
Definition: m_nlte.cc:68
void ArrayOfQuantumIdentifierFromLines(ArrayOfQuantumIdentifier &qid, const ArrayOfArrayOfAbsorptionLines &abs_lines_per_species, const Index &global, const Verbosity &)
WORKSPACE METHOD: ArrayOfQuantumIdentifierFromLines.
Definition: m_nlte.cc:35
void line_irradianceCalcForSingleSpeciesNonOverlappingLinesPseudo2D(Workspace &ws, Matrix &line_irradiance, Tensor3 &line_transmission, const ArrayOfArrayOfSpeciesTag &abs_species, const ArrayOfArrayOfAbsorptionLines &abs_lines_per_species, const EnergyLevelMap &nlte_field, const Tensor4 &vmr_field, const Tensor3 &t_field, const Tensor3 &z_field, const Vector &p_grid, const Vector &refellipsoid, const Tensor3 &surface_props_data, const Agenda &ppath_agenda, const Agenda &iy_main_agenda, const Agenda &iy_space_agenda, const Agenda &iy_surface_agenda, const Agenda &iy_cloudbox_agenda, const Agenda &propmat_clearsky_agenda, const Numeric &df, const Index &nz, const Index &nf, const Numeric &r, const Verbosity &verbosity)
WORKSPACE METHOD: line_irradianceCalcForSingleSpeciesNonOverlappingLinesPseudo2D.
void abs(Sparse &A, const Sparse &B)
Absolute value of sparse matrix elements.
Definition: matpackII.cc:394
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
#define CREATE_OUT2
Definition: messages.h:205
my_basic_string< char > String
The String type for ARTS.
Definition: mystring.h:216
Vector createBij(const ArrayOfArrayOfAbsorptionLines &abs_lines)
Create a Bij object.
Definition: nlte.cc:114
Vector createBji(const Vector &Bij, const ArrayOfArrayOfAbsorptionLines &abs_lines)
Create a Bji object.
Definition: nlte.cc:134
Vector createAij(const ArrayOfArrayOfAbsorptionLines &abs_lines)
Create a Aij object.
Definition: nlte.cc:97
void check_collision_line_identifiers(const ArrayOfQuantumIdentifier &collision_line_identifiers)
Checks that a WSV is OK or throws a run-time error.
Definition: nlte.cc:277
void set_constant_statistical_equilibrium_matrix(MatrixView A, VectorView x, const Numeric &sem_ratio, const Index row)
Set a row of the SEE matrix and level distribution vector to constant.
Definition: nlte.cc:89
void nlte_collision_factorsCalcFromCoeffs(Vector &Cij, Vector &Cji, const ArrayOfArrayOfAbsorptionLines &abs_lines, const ArrayOfArrayOfSpeciesTag &abs_species, const ArrayOfArrayOfGriddedField1 &collision_coefficients, const ArrayOfQuantumIdentifier &collision_line_identifiers, const SpeciesIsotopologueRatios &isotopologue_ratios, const ConstVectorView &vmr, const Numeric &T, const Numeric &P)
Gets collisional factors from coefficients.
Definition: nlte.cc:180
Index find_first_unique_in_lower(const ArrayOfIndex &upper, const ArrayOfIndex &lower) ARTS_NOEXCEPT
Finds a unique lower state if one exists or returns index to last element.
Definition: nlte.cc:268
void nlte_positions_in_statistical_equilibrium_matrix(ArrayOfIndex &upper, ArrayOfIndex &lower, const ArrayOfArrayOfAbsorptionLines &abs_lines, const EnergyLevelMap &nlte_field)
Finds upper and lower states in SEE Matrix.
Definition: nlte.cc:235
void dampened_statistical_equilibrium_equation(MatrixView A, const ConstVectorView &x, const ConstVectorView &Aij, const ConstVectorView &Bij, const ConstVectorView &Bji, const ConstVectorView &Cij, const ConstVectorView &Cji, const ConstVectorView &Jij, const ConstVectorView &Lambda, const ArrayOfIndex &upper, const ArrayOfIndex &lower, const Numeric &total_number_count)
Sets up the solution matrix for linear dampened statistical equilibrium equation.
Definition: nlte.cc:55
void statistical_equilibrium_equation(MatrixView A, const ConstVectorView &Aij, const ConstVectorView &Bij, const ConstVectorView &Bji, const ConstVectorView &Cij, const ConstVectorView &Cji, const ConstVectorView &Jij, const ArrayOfIndex &upper, const ArrayOfIndex &lower)
Sets up the solution matrix for linear statistical equilibrium equation.
Definition: nlte.cc:31
Deep calculations for NLTE.
ArrayOfQuantumIdentifier levels
A logical struct for global quantum numbers with species identifiers.
GlobalState UpperLevel() const
GlobalState LowerLevel() const
This file contains basic functions to handle XML data files.
void xml_read_from_file(const String &filename, T &type, const Verbosity &verbosity)
Reads data from XML file.
Definition: xml_io.h:151