Athena++/Atmosphere
Planetary Atmosphere Simulator
radiation_utils.cpp
Go to the documentation of this file.
1 // C++ header
2 #include <fstream>
3 #include <iostream>
4 #include <sstream>
5 #include <iomanip>
6 
7 // Athena++ header
8 #include "../utils/utils.hpp" // Vectorize
9 #include "../math/core.h" // deg2rad
10 #include "radiation_utils.hpp"
11 
12 /*void Radiation::TotalFlux(AthenaArray<Real>& flux) const
13 {
14  // count how many bands
15  int b = 0;
16  int max_ntau = 0;
17  Radiation const *p = this;
18  while (p != NULL) {
19  max_ntau = std::max(max_ntau, p->ntau);
20  p = p->next;
21  b++;
22  }
23 
24  // check dimension consistancy, reallocate memory if necessary
25  if (flux.GetDim1() != b && flux.GetDim2() < max_ntau) {
26  flux.DeleteAthenaArray();
27  flux.NewAthenaArray(max_ntau, b);
28  }
29 
30  b = 0;
31  p = this;
32  while (p != NULL) {
33  for (int j = 0; j < p->ntau; ++j)
34  flux(j,b) = 0.;
35  for (int i = 0; i < p->nwave; ++i)
36  for (int j = 0; j < p->ntau; ++j)
37  // flup - rfldn - rfldir
38  flux(j,b) += p->weight_[i]*(p->rad(i,j,2) - p->rad(i,j,1) - p->rad(i,j,0));
39  p = p->next;
40  b++;
41  }
42 }
43 
44 void Radiation::SumTotalFlux(AthenaArray<Real>& tflux, AthenaArray<Real>& disort_flux, int gk, int gj) const
45 {
46  // count how many bands
47  int b = 0;
48  int max_ntau = 0;
49  Radiation const *p = this;
50  while (p != NULL) {
51  max_ntau = std::max(max_ntau, p->ntau);
52  p = p->next;
53  b++;
54  }
55  AthenaArray<Real> flux;
56  flux.NewAthenaArray(max_ntau, b);
57  b = 0;
58  p = this;
59  while (p != NULL) {
60  for (int j = 0; j < p->ntau; ++j)
61  flux(j,b) = 0.;
62  for (int j = 0; j < p->ntau+2*(NGHOST); ++j){
63  disort_flux(0,gk,gj,j) = 0.;
64  disort_flux(1,gk,gj,j) = 0.;
65  disort_flux(2,gk,gj,j) = 0.;
66  }
67  for (int i = 0; i < p->nwave; ++i)
68  for (int j = 0; j < p->ntau; ++j){
69  // flup - rfldn - rfldir
70  disort_flux(0,gk,gj,j+NGHOST) += p->weight_[i]*p->rad(i,j,2);
71  disort_flux(1,gk,gj,j+NGHOST) += p->weight_[i]*p->rad(i,j,1);
72  disort_flux(2,gk,gj,j+NGHOST) += p->weight_[i]*p->rad(i,j,0);
73  flux(j,b) += p->weight_[i]*(p->rad(i,j,2) - p->rad(i,j,1) - p->rad(i,j,0));
74  }
75  p = p->next;
76  b++;
77  }
78 
79  // sum; thl
80  for (int j =0; j < max_ntau; ++j){
81  tflux(j) = 0.;
82  }
83  for (int j =0; j <max_ntau; ++j)
84  for (int bb =0; bb <b; ++bb){
85  tflux(j) += flux(max_ntau-j-1,bb);
86  }
87 }
88 
89 void Radiation::WriteTopFlux(std::string fname) const
90 {
91  std::cout << "Top flux written into file: " << fname << std::endl;
92  std::ofstream out(fname.c_str(), std::ios::out);
93  out << std::left << std::setw(20) << std::setprecision(8) << "# Wavenumber(cm-1)"
94  << std::left << std::setw(20) << "Top rfldir[]"
95  << std::left << std::setw(20) << "Top rfldn[]"
96  << std::left << std::setw(20) << "Top flup[]"
97  << std::endl;
98 
99  Radiation const *p = this;
100  while (p != NULL) {
101  for (int i = 0; i < p->nwave; ++i)
102  out << std::left << std::setw(20) << p->wave[i]
103  << std::left << std::setw(20) << p->rad(i,0,0)
104  << std::left << std::setw(20) << p->rad(i,0,1)
105  << std::left << std::setw(20) << p->rad(i,0,2)
106  << std::endl;
107  p = p->next;
108  }
109 }
110 
111 void Radiation::WriteOpticalDepth(std::string fname) const
112 {
113  std::cout << "Optical depth written into file: " << fname << std::endl;
114  std::ofstream out(fname.c_str(), std::ios::out);
115 
116  Radiation const *p = this;
117  while (p != NULL) {
118  out << std::left << "# Band: " << p->myname << std::endl;
119  out << std::setw(40) << "# Number of wavenumbers : " << p->nwave << std::endl;
120  out << std::setw(40) << "# Number of levels: " << p->nlevel << std::endl;
121  out << std::right;
122  for (int i = 0; i < p->nlevel; ++i)
123  out << std::setw(16) << std::setprecision(8) << p->level[i]/1.E3;
124  out << std::endl;
125 
126  for (int i = 0; i < p->nwave; ++i) {
127  out << std::setw(16) << std::setprecision(8) << p->wave[i];
128  for (int j = 0; j < p->nlevel - 1; ++j)
129  out << std::setw(16) << std::setprecision(8) << p->tau(ITAU,i,j);
130  out << std::endl;
131  }
132  p = p->next;
133  }
134 }
135 
136 void Radiation::WriteTopRadiance(std::string fname) const
137 {
138  std::cout << "Top radiance written into file: " << fname << std::endl;
139  std::ofstream out(fname.c_str(), std::ios::out);
140  out << std::left << std::setw(20) << std::setprecision(8)
141  << "# Wavenumber/Frequency" << std::endl;
142 
143  Radiation const *p = this;
144  while (p != NULL) {
145  for (int i = 0; i < p->nwave; ++i) {
146  out << std::left << std::setw(20) << p->wave[i];
147  for (int j = 0; j < p->nphi; ++j)
148  for (int k = 0; k < p->numu; ++k)
149  out << std::left << std::setw(20) << p->uu(i,j,0,k);
150  out << std::endl;
151  }
152  p = p->next;
153  }
154 }
155 
156 void WriteHeatingRate(std::string fname, AthenaArray<Real> const& flux,
157  AthenaArray<Real> const& hrate, Real const* level)
158 {
159  int nband = flux.GetDim1();
160  int nlevel = flux.GetDim2();
161 
162  // print mean intensity, flux divergence and heating rate
163  std::cout << "Heating rate written into file: " << fname << std::endl;
164  std::ofstream out(fname.c_str(), std::ios::out);
165  out << std::left << std::setw(6) << "Level"
166  << std::left << std::setw(12) << "Height [km]";
167 
168  for (int b = 0; b < nband; ++b) {
169  char bname[80];
170  sprintf(bname, "B%d flux [w/m^2]", b + 1);
171  out << std::left << std::setw(16) << bname
172  << std::left << std::setw(20) << "Heating rate [w/kg]";
173  }
174 
175  if (nband > 1) {
176  out << std::left << std::setw(16) << "Flux [w/m^2]"
177  << std::left << std::setw(20) << "Heating rate [w/kg]"
178  << std::endl;
179  } else {
180  out << std::endl;
181  }
182 
183  for (int i = 0; i < nlevel; ++i) {
184  out << std::left << std::setw(6) << i + 1
185  << std::left << std::setw(12) << level[i]/1.E3;
186 
187  Real total_flux = 0., total_hrate = 0.;
188  for (int b = 0; b < nband; ++b) {
189  out << std::left << std::setw(16) << std::setprecision(8) << flux(i,b)
190  << std::left << std::setw(20) << std::setprecision(8) << hrate(i,b);
191  total_flux += flux(i,b);
192  total_hrate += hrate(i,b);
193  }
194 
195  if (nband > 1) {
196  out << std::left << std::setw(16) << std::setprecision(8) << total_flux
197  << std::left << std::setw(20) << std::setprecision(8) << total_hrate
198  << std::endl;
199  } else {
200  out << std::endl;
201  }
202  }
203 }*/
204 
205 void setRadiationFlags(uint64_t *flags, std::string str)
206 {
207  std::stringstream msg;
208  std::vector<std::string> dstr = Vectorize<std::string>(str.c_str(), " ,");
209  for (int i = 0; i < dstr.size(); ++i) {
210  if (dstr[i] == "static") {
211  *flags &= !RadiationFlags::Dynamic;
212  } else if (dstr[i] == "dynamic") {
213  *flags |= RadiationFlags::Dynamic;
214  } else if (dstr[i] == "bin") {
215  *flags &= !RadiationFlags::LineByLine;
216  } else if (dstr[i] == "lbl") {
217  *flags |= RadiationFlags::LineByLine;
218  } else if (dstr[i] == "ck") {
219  *flags |= RadiationFlags::CorrelatedK;
220  } else if (dstr[i] == "planck") {
221  *flags |= RadiationFlags::Planck;
222  } else if (dstr[i] == "star") {
223  *flags |= RadiationFlags::Star;
224  } else if (dstr[i] == "spher") {
225  *flags |= RadiationFlags::Sphere;
226  } else if (dstr[i] == "only") {
227  *flags |= RadiationFlags::FluxOnly;
228  } else {
229  msg << "### FATAL ERROR in function setRadiatino"
230  << std::endl << "flag:" << dstr[i] << "unrecognized"
231  << std::endl;
232  ATHENA_ERROR(msg);
233  }
234 
235  // check flags consistency
236  if ((*flags & RadiationFlags::LineByLine) && (*flags & RadiationFlags::CorrelatedK)) {
237  msg << "### FATAL ERROR in function setRadiation"
238  << std::endl << "ck cannot be used with lbl."
239  << std::endl;
240  }
241  }
242 }
243 
244 void getPhaseHenyeyGreenstein(Real *pmom, int iphas, Real gg, int npmom)
245 {
246  pmom[0] = 1.;
247  for (int k =1; k <npmom; k++)
248  pmom[k] = 0.;
249 
250  switch(iphas){
251  case 0: // HENYEY_GREENSTEIN
252  if (gg <= -1. || gg >= 1.)
253  std::cout<<"getmom--bad input variable gg"<<std::endl;
254  for (int k = 1; k <= npmom; k++)
255  pmom[k] = pow(gg,(Real)k);
256  break;
257 
258  case 1: // RAYLEIGH
259  pmom[2] = 0.1;
260  break;
261  }
262 }
263 
264 void packSpectralProperties(Real *buf, Real const *tau, Real const *ssa, Real const* pmom, int nlayer, int npmom)
265 {
266  for (int i = 0; i < nlayer; ++i)
267  *(buf++) = tau[i];
268  for (int i = 0; i < nlayer; ++i)
269  *(buf++) = ssa[i];
270  for (int i = 0; i < nlayer; ++i)
271  for (int j = 0; j < npmom; ++j)
272  *(buf++) = pmom[i*npmom+j];
273 }
274 
275 void unpackSpectralProperties(Real *tau, Real *ssa, Real *pmom, Real const *buf, int slyr, int npmom, int nblocks, int npmom_max)
276 {
277  npmom_max = std::max(npmom_max, npmom);
278  for (int n = 0; n < nblocks; ++n) {
279  for (int i = 0; i < slyr; ++i)
280  tau[n*slyr + i] = *(buf++);
281  for (int i = 0; i < slyr; ++i)
282  ssa[n*slyr + i] = *(buf++);
283  for (int i = 0; i < slyr; ++i) {
284  for (int j = 0; j < npmom; ++j)
285  pmom[n*slyr*npmom_max + i*npmom_max + j] = *(buf++);
286  for (int j = npmom; j < npmom_max; ++j)
287  pmom[n*slyr*npmom_max + i*npmom_max + j] = 0.;
288  }
289  }
290 }
291 
292 void readRadiationDirections(std::vector<Direction> &ray, std::string str)
293 {
294  std::vector<std::string> dstr = Vectorize<std::string>(str.c_str());
295  int nray = dstr.size();
296  ray.resize(nray);
297 
298  for (int i = 0; i < nray; ++i) {
299  ray[i].phi = 0.;
300  sscanf(dstr[i].c_str(), "(%lf,%lf)", &ray[i].mu, &ray[i].phi);
301  ray[i].mu = cos(deg2rad(ray[i].mu));
302  ray[i].phi = deg2rad(ray[i].phi);
303  }
304 }
double Real
Definition: athena.hpp:29
double deg2rad(double phi)
Definition: core.h:43
double max(double x1, double x2, double x3)
Definition: core.h:18
const uint64_t Sphere
Definition: radiation.hpp:39
const uint64_t Dynamic
Definition: radiation.hpp:34
const uint64_t Star
Definition: radiation.hpp:38
const uint64_t LineByLine
Definition: radiation.hpp:35
const uint64_t FluxOnly
Definition: radiation.hpp:40
const uint64_t Planck
Definition: radiation.hpp:37
const uint64_t CorrelatedK
Definition: radiation.hpp:36
void getPhaseHenyeyGreenstein(Real *pmom, int iphas, Real gg, int npmom)
void readRadiationDirections(std::vector< Direction > &ray, std::string str)
void packSpectralProperties(Real *buf, Real const *tau, Real const *ssa, Real const *pmom, int nlayer, int npmom)
void unpackSpectralProperties(Real *tau, Real *ssa, Real *pmom, Real const *buf, int slyr, int npmom, int nblocks, int npmom_max)
void setRadiationFlags(uint64_t *flags, std::string str)