ARTS 2.5.11 (git: 725533f0)
m_nlte.cc
Go to the documentation of this file.
1
9#include "absorption.h"
10#include "arts.h"
11#include "auto_md.h"
12#include "lin_alg.h"
13#include "nlte.h"
14#include "xml_io.h"
15
16/* Workspace method: Doxygen documentation will be auto-generated */
19 const ArrayOfArrayOfAbsorptionLines& abs_lines_per_species,
20 const Index& global,
21 const Verbosity&) {
22 // Defined only as output not input so resizing
23 qid.resize(0);
24
25 QuantumIdentifier lower, upper;
26
27 // For all lines' upper and lower energy levels
28 for (const auto& lines: abs_lines_per_species) {
29 for (const auto& band: lines) {
30 for (Index k=0; k<band.NumLines() and (global ? (k==0) : false); k++) {
31 if (global) {
32 lower = band.quantumidentity.LowerLevel();
33 upper = band.quantumidentity.UpperLevel();
34 } else {
35 auto x = band.QuantumIdentityOfLine(k);
36 lower = x.LowerLevel();
37 upper = x.UpperLevel();
38 }
39
40 if (std::none_of(qid.begin(), qid.end(), [&](auto& x){return x == lower;}))
41 qid.push_back(lower);
42 if (std::none_of(qid.begin(), qid.end(), [&](auto& x){return x == upper;}))
43 qid.push_back(upper);
44 }
45 }
46 }
47}
48
49/* Workspace method: Doxygen documentation will be auto-generated */
51 const Numeric& scale,
52 const Verbosity&) {
53 nlte_field.value *= scale;
54}
55
56/* Workspace method: Doxygen documentation will be auto-generated */
58 Workspace& ws,
59 EnergyLevelMap& nlte_field,
60 const ArrayOfArrayOfSpeciesTag& abs_species,
61 const ArrayOfArrayOfAbsorptionLines& abs_lines_per_species,
62 const ArrayOfArrayOfGriddedField1& collision_coefficients,
63 const ArrayOfQuantumIdentifier& collision_line_identifiers,
64 const SpeciesIsotopologueRatios& isotopologue_ratios,
65 const Agenda& iy_main_agenda,
66 const Agenda& ppath_agenda,
67 const Agenda& iy_space_agenda,
68 const Agenda& iy_surface_agenda,
69 const Agenda& iy_cloudbox_agenda,
70 const Agenda& propmat_clearsky_agenda,
71 const Agenda& /*water_p_eq_agenda*/,
72 const Tensor4& vmr_field,
73 const Tensor3& t_field,
74 const Tensor3& z_field,
75 const Vector& p_grid,
76 const Index& atmosphere_dim,
77 const Vector& refellipsoid,
78 const Tensor3& surface_props_data,
79 const Index& nlte_do,
80 const Numeric& df,
81 const Numeric& convergence_limit,
82 const Index& nz,
83 const Index& nf,
84 const Index& dampened,
85 const Index& iteration_limit,
86 const Verbosity& verbosity)
87{
89
90 ARTS_USER_ERROR_IF (not nlte_do, "Must be set to do NLTE");
91 ARTS_USER_ERROR_IF (nlte_field.value.empty(),
92 "Error in NLTE field, it is empty");
93
94 Matrix line_irradiance;
95 Tensor3 line_transmission;
96
97 const Index nlevels = nlte_field.levels.nelem(), np = p_grid.nelem();
98 ARTS_USER_ERROR_IF (nlevels < 5,
99 "Must have more than a four levels");
100
101 ARTS_USER_ERROR_IF (atmosphere_dim not_eq 1,
102 "Only for 1D atmosphere");
103
104 const Index nlines = nelem(abs_lines_per_species);
105 ARTS_USER_ERROR_IF (nlevels >= nlines,
106 "Bad number of lines... overlapping lines in nlte_level_identifiers?");
107
108 // Create basic compute vectors
109 const Vector Aij = createAij(abs_lines_per_species);
110 const Vector Bij = createBij(abs_lines_per_species);
111 const Vector Bji = createBji(Bij, abs_lines_per_species);
112 Vector Cij(nlines), Cji(nlines);
113
114 ArrayOfIndex upper, lower;
116 upper, lower, abs_lines_per_species, nlte_field);
117 const Index unique = find_first_unique_in_lower(upper, lower);
118
119 // Compute arrays
120 Matrix SEE(nlevels, nlevels, 0.0);
121 Vector r(nlevels, 0.0), x(nlevels, 0.0);
122 Numeric max_change = convergence_limit + 1;
123
124 Index i = 0;
125 while (i < iteration_limit and max_change > convergence_limit) {
126 // Reset change
127 max_change = 0.0;
128
129 // Compute radiation and transmission
131 ws,
132 line_irradiance,
133 line_transmission,
134 abs_species,
135 abs_lines_per_species,
136 nlte_field,
137 vmr_field,
138 t_field,
139 z_field,
140 p_grid,
141 refellipsoid,
142 surface_props_data,
143 ppath_agenda,
144 iy_main_agenda,
145 iy_space_agenda,
146 iy_surface_agenda,
147 iy_cloudbox_agenda,
148 propmat_clearsky_agenda,
149 df,
150 nz,
151 nf,
152 1.0,
153 verbosity);
154
155 for (Index ip = 0; ip < np; ip++) {
156 r = nlte_field.value(joker, ip, 0, 0);
158 Cji,
159 abs_lines_per_species,
160 abs_species,
161 collision_coefficients,
162 collision_line_identifiers,
163 isotopologue_ratios,
164 vmr_field(joker, ip, 0, 0),
165 t_field(ip, 0, 0),
166 p_grid[ip]);
167
168
169 if (dampened)
171 SEE,
172 r,
173 Aij,
174 Bij,
175 Bji,
176 Cij,
177 Cji,
178 line_irradiance(joker, ip),
179 line_transmission(0, joker, ip),
180 upper,
181 lower);
182 else
184 Aij,
185 Bij,
186 Bji,
187 Cij,
188 Cji,
189 line_irradiance(joker, ip),
190 upper,
191 lower);
192
193 set_constant_statistical_equilibrium_matrix(SEE, x, sum(r), unique);
194 solve(nlte_field.value(joker, ip, 0, 0), SEE, x);
195
196 for (Index il = 0; il < nlevels; il++) {
197 max_change =
198 max(abs(nlte_field.value(il, ip, 0, 0) - r[il]) / r[il], max_change);
199 }
200 }
201 i++;
202 }
203
204 if (i < iteration_limit)
205 out2 << "Converged NLTE ratios (within convergence_limit) returned after "
206 << i << " iterations\n";
207 else
208 out2
209 << "No convergence of NLTE ratios (within convergence_limit) returned even after "
210 << iteration_limit << " iterations\n";
211}
212
213/* Workspace method: Doxygen documentation will be auto-generated */
215 ArrayOfArrayOfGriddedField1& collision_coefficients,
216 ArrayOfQuantumIdentifier& collision_line_identifiers,
217 const ArrayOfArrayOfSpeciesTag& abs_species,
218 const String& basename,
219 const Verbosity& verbosity) {
220 // Standard ARTS basename-assumption
221 String tmp_basename = basename, filename;
222 if (basename.length() && basename[basename.length() - 1] != '/')
223 tmp_basename += ".";
224
225 // Read the identification file
226 filename = tmp_basename + "qid.xml";
227 xml_read_from_file(filename, collision_line_identifiers, verbosity);
228 check_collision_line_identifiers(collision_line_identifiers);
229
230 // Inner array size has to be this constantly
231 const Index n = collision_line_identifiers.nelem();
232
233 // Set species dimensions and fill the array
234 collision_coefficients.resize(abs_species.nelem());
235 for (Index i = 0; i < collision_coefficients.nelem(); i++) {
237
238 // Read the file for a species and check that the size is correct of the array
239 filename = tmp_basename + String(Species::toShortName(abs_species[i].Species())) + ".xml";
240 xml_read_from_file(filename, aogf1, verbosity);
241 ARTS_USER_ERROR_IF (aogf1.nelem() not_eq n,
242 "Mismatch between collision_line_identifiers and some collision_coefficients");
243 collision_coefficients[i] = aogf1;
244 }
245}
246
247/* Workspace method: Doxygen documentation will be auto-generated */
248void nlteOff(Index& nlte_do,
249 EnergyLevelMap& nlte_field,
250 ArrayOfQuantumIdentifier& nlte_level_identifiers,
251 const Verbosity&) {
252 nlte_do = 0;
253 nlte_field = EnergyLevelMap();
254 nlte_level_identifiers.resize(0);
255}
Declarations required for the calculation of absorption coefficients.
base max(const Array< base > &x)
Max function.
Definition array.h:128
The global header file for ARTS.
The Agenda class.
Index nelem() const ARTS_NOEXCEPT
Definition array.h:75
Workspace class.
void solve(VectorView w, const CovarianceMatrix &A, ConstVectorView v)
#define ARTS_USER_ERROR_IF(condition,...)
Definition debug.h:137
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:57
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:214
void nlteOff(Index &nlte_do, EnergyLevelMap &nlte_field, ArrayOfQuantumIdentifier &nlte_level_identifiers, const Verbosity &)
WORKSPACE METHOD: nlteOff.
Definition m_nlte.cc:248
void nlte_fieldRescalePopulationLevels(EnergyLevelMap &nlte_field, const Numeric &scale, const Verbosity &)
WORKSPACE METHOD: nlte_fieldRescalePopulationLevels.
Definition m_nlte.cc:50
void ArrayOfQuantumIdentifierFromLines(ArrayOfQuantumIdentifier &qid, const ArrayOfArrayOfAbsorptionLines &abs_lines_per_species, const Index &global, const Verbosity &)
WORKSPACE METHOD: ArrayOfQuantumIdentifierFromLines.
Definition m_nlte.cc:17
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.
#define CREATE_OUT2
Definition messages.h:188
my_basic_string< char > String
The String type for ARTS.
Definition mystring.h:199
Vector createBij(const ArrayOfArrayOfAbsorptionLines &abs_lines)
Create a Bij object.
Definition nlte.cc:96
Vector createBji(const Vector &Bij, const ArrayOfArrayOfAbsorptionLines &abs_lines)
Create a Bji object.
Definition nlte.cc:116
Vector createAij(const ArrayOfArrayOfAbsorptionLines &abs_lines)
Create a Aij object.
Definition nlte.cc:79
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:259
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:71
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:162
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:250
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:217
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:37
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:13
Deep calculations for NLTE.
ArrayOfQuantumIdentifier levels
A logical struct for global quantum numbers with species identifiers.
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:134