ARTS  2.4.0(git:4fb77825)
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 
33 /* Workspace method: Doxygen documentation will be auto-generated */
37  const Index& global,
38  const Verbosity&) {
39  // Defined only as output not input so resizing
40  qid.resize(0);
41 
42  QuantumIdentifier lower, upper;
43 
44  // For all lines' upper and lower energy levels
45  for (const auto& lines: abs_lines_per_species) {
46  for (const auto& band: lines) {
47  for (Index k=0; k<band.NumLines() and (global ? (k==0) : false); k++) {
48  if (global) {
49  lower = band.QuantumIdentity().LowerQuantumId();
50  upper = band.QuantumIdentity().UpperQuantumId();
51  } else {
52  auto x = band.QuantumIdentityOfLine(k);
53  lower = x.LowerQuantumId();
54  upper = x.UpperQuantumId();
55  }
56 
57  const bool canbeinlower = lower.any_quantumnumbers(),
58  canbeinupper = upper.any_quantumnumbers();
59 
60  // Test if the level has already been treated
61  const Index n = qid.nelem();
62  bool inlower = false, inupper = false;
63  for (Index i = 0; i < n; i++) {
64  if (not inlower and canbeinlower)
65  if (qid[i] == lower) inlower = true;
66  if (not inupper and canbeinupper)
67  if (qid[i] == upper) inupper = true;
68  }
69 
70  // If the level has not been treated and has any quantum numbers, then store it
71  if (not inlower and canbeinlower) qid.push_back(lower);
72  if (not inupper and canbeinupper)
73  if (not(lower == upper)) qid.push_back(upper);
74  }
75  }
76  }
77 }
78 
79 /* Workspace method: Doxygen documentation will be auto-generated */
81  const Numeric& scale,
82  const Verbosity&) {
83  nlte_field.Data() *= scale;
84 }
85 
86 /* Workspace method: Doxygen documentation will be auto-generated */
88  Workspace& ws,
95  const Agenda& iy_main_agenda,
96  const Agenda& ppath_agenda,
97  const Agenda& iy_space_agenda,
101  const Agenda& /*water_p_eq_agenda*/,
102  const Tensor4& vmr_field,
103  const Tensor3& t_field,
104  const Tensor3& z_field,
105  const Vector& p_grid,
106  const Index& atmosphere_dim,
107  const Vector& refellipsoid,
109  const Index& nlte_do,
110  const Numeric& df,
111  const Numeric& convergence_limit,
112  const Index& nz,
113  const Index& nf,
114  const Index& dampened,
115  const Index& iteration_limit,
116  const Verbosity& verbosity)
117 {
118  CREATE_OUT2;
119 
120  if (not nlte_do) throw std::runtime_error("Must be set to do NLTE");
121  if (nlte_field.Data().empty())
122  throw std::runtime_error("Error in NLTE field, it is empty");
123 
126 
127  const Index nlevels = nlte_field.Levels().nelem(), np = p_grid.nelem();
128  if (nlevels < 5)
129  throw std::runtime_error("Must have more than a four levels");
130 
131  if (atmosphere_dim not_eq 1)
132  throw std::runtime_error("Only for 1D atmosphere");
133 
135  if (nlevels >= nlines)
136  throw std::runtime_error(
137  "Bad number of lines... overlapping lines in nlte_level_identifiers?");
138 
139  // Create basic compute vectors
142  const Vector Bji = createBji(Bij, abs_lines_per_species);
143  Vector Cij(nlines), Cji(nlines);
144 
145  ArrayOfIndex upper, lower;
147  upper, lower, abs_lines_per_species, nlte_field);
148  const Index unique = find_first_unique_in_lower(upper, lower);
149 
150  // Compute arrays
151  Matrix SEE(nlevels, nlevels, 0.0);
152  Vector r(nlevels, 0.0), x(nlevels, 0.0);
153  Numeric max_change = convergence_limit + 1;
154 
155  Index i = 0;
156  while (i < iteration_limit and max_change > convergence_limit) {
157  // Reset change
158  max_change = 0.0;
159 
160  // //Compute radiation and transmission
162  ws,
165  abs_species,
167  nlte_field,
168  vmr_field,
169  t_field,
170  z_field,
171  p_grid,
172  refellipsoid,
174  ppath_agenda,
180  df,
181  nz,
182  nf,
183  1.0,
184  verbosity);
185 
186  for (Index ip = 0; ip < np; ip++) {
187  r = nlte_field.Data()(joker, ip, 0, 0);
189  Cji,
191  abs_species,
195  vmr_field(joker, ip, 0, 0),
196  t_field(ip, 0, 0),
197  p_grid[ip]);
198 
199 
200  if (dampened)
202  SEE,
203  r,
204  Aij,
205  Bij,
206  Bji,
207  Cij,
208  Cji,
209  line_irradiance(joker, ip),
210  line_transmission(0, joker, ip),
211  upper,
212  lower);
213  else
215  Aij,
216  Bij,
217  Bji,
218  Cij,
219  Cji,
220  line_irradiance(joker, ip),
221  upper,
222  lower);
223 
224  set_constant_statistical_equilibrium_matrix(SEE, x, r.sum(), unique);
225  solve(nlte_field.Data()(joker, ip, 0, 0), SEE, x);
226 
227  for (Index il = 0; il < nlevels; il++) {
228  max_change =
229  max(abs(nlte_field.Data()(il, ip, 0, 0) - r[il]) / r[il], max_change);
230  }
231  }
232  i++;
233  }
234 
235  if (i < iteration_limit)
236  out2 << "Converged NLTE ratios (within convergence_limit) returned after "
237  << i << " iterations\n";
238  else
239  out2
240  << "No convergence of NLTE ratios (within convergence_limit) returned even after "
241  << iteration_limit << " iterations\n";
242 }
243 
244 /* Workspace method: Doxygen documentation will be auto-generated */
249  const String& basename,
250  const Verbosity& verbosity) {
251  // Standard ARTS basename-assumption
252  String tmp_basename = basename, filename;
253  if (basename.length() && basename[basename.length() - 1] != '/')
254  tmp_basename += ".";
255 
256  // Read the identification file
257  filename = tmp_basename + "qid.xml";
260 
261  // Inner array size has to be this constantly
262  const Index n = collision_line_identifiers.nelem();
263 
264  // Set species dimensions and fill the array
265  collision_coefficients.resize(abs_species.nelem());
266  for (Index i = 0; i < collision_coefficients.nelem(); i++) {
267  ArrayOfGriddedField1 aogf1;
268 
269  // Read the file for a species and check that the size is correct of the array
270  filename = tmp_basename + abs_species[i][0].SpeciesNameMain() + ".xml";
271  xml_read_from_file(filename, aogf1, verbosity);
272  if (aogf1.nelem() not_eq n)
273  throw std::runtime_error(
274  "Mismatch between collision_line_identifiers and some collision_coefficients");
275  collision_coefficients[i] = aogf1;
276  }
277 }
278 
279 /* Workspace method: Doxygen documentation will be auto-generated */
283  const Verbosity&) {
284  nlte_do = 0;
286  nlte_level_identifiers.resize(0);
287 }
ARTS::Var::collision_coefficients
ArrayOfArrayOfGriddedField1 collision_coefficients(Workspace &ws) noexcept
Definition: autoarts.h:2789
Matrix
The Matrix class.
Definition: matpackI.h:1193
ARTS::Var::atmosphere_dim
Index atmosphere_dim(Workspace &ws) noexcept
Definition: autoarts.h:2510
check_collision_line_identifiers
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:286
ARTS::Var::z_field
Tensor3 z_field(Workspace &ws) noexcept
Definition: autoarts.h:7690
nlte_collision_factorsCalcFromCoeffs
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 SpeciesAuxData &isotopologue_ratios, const ConstVectorView vmr, const Numeric &T, const Numeric &P)
Gets collisional factors from coefficients.
Definition: nlte.cc:182
QuantumIdentifier
Class to identify and match lines by their quantum numbers.
Definition: quantum.h:390
auto_md.h
ARTS::Var::iy_cloudbox_agenda
Agenda iy_cloudbox_agenda(Workspace &ws) noexcept
Definition: autoarts.h:3742
absorption.h
Declarations required for the calculation of absorption coefficients.
QuantumIdentifier::UpperQuantumId
constexpr QuantumIdentifier UpperQuantumId() const noexcept
Return a quantum identifer as if it wants to match to upper energy level.
Definition: quantum.h:577
Tensor3
The Tensor3 class.
Definition: matpackIII.h:339
QuantumIdentifier::any_quantumnumbers
bool any_quantumnumbers() const
Check if there are any quantum numbers defined.
Definition: quantum.cc:392
find_first_unique_in_lower
Index find_first_unique_in_lower(const ArrayOfIndex &upper, const ArrayOfIndex &lower) noexcept
Finds a unique lower state if one exists or returns index to last element.
Definition: nlte.cc:277
joker
const Joker joker
abs
#define abs(x)
Definition: legacy_continua.cc:20626
ARTS::Var::verbosity
Verbosity verbosity(Workspace &ws) noexcept
Definition: autoarts.h:7112
ARTS::Var::nlte_level_identifiers
ArrayOfQuantumIdentifier nlte_level_identifiers(Workspace &ws) noexcept
Definition: autoarts.h:4639
ARTS::Var::abs_lines_per_species
ArrayOfArrayOfAbsorptionLines abs_lines_per_species(Workspace &ws) noexcept
Definition: autoarts.h:2022
ARTS::Var::line_irradiance
Matrix line_irradiance(Workspace &ws) noexcept
Definition: autoarts.h:4005
nlteOff
void nlteOff(Index &nlte_do, EnergyLevelMap &nlte_field, ArrayOfQuantumIdentifier &nlte_level_identifiers, const Verbosity &)
WORKSPACE METHOD: nlteOff.
Definition: m_nlte.cc:280
ARTS::Var::collision_line_identifiers
ArrayOfQuantumIdentifier collision_line_identifiers(Workspace &ws) noexcept
Definition: autoarts.h:2797
CREATE_OUT2
#define CREATE_OUT2
Definition: messages.h:206
Tensor4
The Tensor4 class.
Definition: matpackIV.h:421
ARTS::Var::iy_space_agenda
Agenda iy_space_agenda(Workspace &ws) noexcept
Definition: autoarts.h:3801
SpeciesAuxData
Auxiliary data for isotopologues.
Definition: absorption.h:217
Agenda
The Agenda class.
Definition: agenda_class.h:44
createBji
Vector createBji(const Vector &Bij, const ArrayOfArrayOfAbsorptionLines &abs_lines)
Create a Bji object.
Definition: nlte.cc:135
createBij
Vector createBij(const ArrayOfArrayOfAbsorptionLines &abs_lines)
Create a Bij object.
Definition: nlte.cc:113
ARTS::Var::nlte_field
EnergyLevelMap nlte_field(Workspace &ws) noexcept
Definition: autoarts.h:4604
Array< QuantumIdentifier >
line_irradianceCalcForSingleSpeciesNonOverlappingLinesPseudo2D
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.
Definition: m_radiation_field.cc:37
Absorption::nelem
Index nelem(const Lines &l)
Number of lines.
Definition: absorptionlines.h:1820
ArrayOfQuantumIdentifierFromLines
void ArrayOfQuantumIdentifierFromLines(ArrayOfQuantumIdentifier &qid, const ArrayOfArrayOfAbsorptionLines &abs_lines_per_species, const Index &global, const Verbosity &)
WORKSPACE METHOD: ArrayOfQuantumIdentifierFromLines.
Definition: m_nlte.cc:34
xml_read_from_file
void xml_read_from_file(const String &filename, T &type, const Verbosity &verbosity)
Reads data from XML file.
Definition: xml_io.cc:901
createAij
Vector createAij(const ArrayOfArrayOfAbsorptionLines &abs_lines)
Create a Aij object.
Definition: nlte.cc:96
my_basic_string< char >
ARTS::Var::vmr_field
Tensor4 vmr_field(Workspace &ws) noexcept
Definition: autoarts.h:7130
ARTS::Var::abs_species
ArrayOfArrayOfSpeciesTag abs_species(Workspace &ws) noexcept
Definition: autoarts.h:2157
ARTS::Var::p_grid
Vector p_grid(Workspace &ws) noexcept
Definition: autoarts.h:4763
ARTS::Var::surface_props_data
Tensor3 surface_props_data(Workspace &ws) noexcept
Definition: autoarts.h:6745
Numeric
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
lin_alg.h
Linear algebra functions.
Verbosity
Definition: messages.h:49
QuantumIdentifier::LowerQuantumId
constexpr QuantumIdentifier LowerQuantumId() const noexcept
Return a quantum identifer as if it wants to match to lower energy level.
Definition: quantum.h:582
ARTS::Var::refellipsoid
Vector refellipsoid(Workspace &ws) noexcept
Definition: autoarts.h:5519
EnergyLevelMap
Definition: energylevelmap.h:60
max
#define max(a, b)
Definition: legacy_continua.cc:20629
ARTS::Var::line_transmission
Tensor3 line_transmission(Workspace &ws) noexcept
Definition: autoarts.h:4014
ARTS::Var::iy_surface_agenda
Agenda iy_surface_agenda(Workspace &ws) noexcept
Definition: autoarts.h:3808
ARTS::Var::nlte_do
Index nlte_do(Workspace &ws) noexcept
Definition: autoarts.h:4567
nlines
#define nlines
Definition: legacy_continua.cc:22187
ARTS::Var::ppath_agenda
Agenda ppath_agenda(Workspace &ws) noexcept
Definition: autoarts.h:5146
nlte_positions_in_statistical_equilibrium_matrix
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:244
ARTS::Var::iy_main_agenda
Agenda iy_main_agenda(Workspace &ws) noexcept
Definition: autoarts.h:3794
nlte.h
Deep calculations for NLTE.
set_constant_statistical_equilibrium_matrix
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:88
ARTS::Var::t_field
Tensor3 t_field(Workspace &ws) noexcept
Definition: autoarts.h:6947
Workspace
Workspace class.
Definition: workspace_ng.h:40
statistical_equilibrium_equation
void statistical_equilibrium_equation(MatrixView A, ConstVectorView Aij, ConstVectorView Bij, ConstVectorView Bji, ConstVectorView Cij, ConstVectorView Cji, ConstVectorView Jij, const ArrayOfIndex &upper, const ArrayOfIndex &lower)
Sets up the solution matrix for linear statistical equilibrium equation.
Definition: nlte.cc:30
ARTS::Var::isotopologue_ratios
SpeciesAuxData isotopologue_ratios(Workspace &ws) noexcept
Definition: autoarts.h:3669
ARTS::Var::propmat_clearsky_agenda
Agenda propmat_clearsky_agenda(Workspace &ws) noexcept
Definition: autoarts.h:5405
ARTS::Group::EnergyLevelMap
EnergyLevelMap EnergyLevelMap
Definition: autoarts.h:78
ARTS::Var::x
Vector x(Workspace &ws) noexcept
Definition: autoarts.h:7346
Index
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
Vector
The Vector class.
Definition: matpackI.h:860
Array::nelem
Index nelem() const
Number of elements.
Definition: array.h:195
nlte_fieldForSingleSpeciesNonOverlappingLines
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 SpeciesAuxData &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:87
collision_coefficientsFromSplitFiles
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:245
nlte_fieldRescalePopulationLevels
void nlte_fieldRescalePopulationLevels(EnergyLevelMap &nlte_field, const Numeric &scale, const Verbosity &)
WORKSPACE METHOD: nlte_fieldRescalePopulationLevels.
Definition: m_nlte.cc:80
arts.h
The global header file for ARTS.
dampened_statistical_equilibrium_equation
void dampened_statistical_equilibrium_equation(MatrixView A, ConstVectorView x, ConstVectorView Aij, ConstVectorView Bij, ConstVectorView Bji, ConstVectorView Cij, ConstVectorView Cji, ConstVectorView Jij, 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:54
solve
void solve(VectorView w, const CovarianceMatrix &A, ConstVectorView v)
Definition: covariance_matrix.cc:594