Athena++/Atmosphere
Planetary Atmosphere Simulator
radiation.cpp
Go to the documentation of this file.
1 // C/C++ headers
2 #include <sstream>
3 #include <stdexcept>
4 
5 // Athena++ header
6 //#include "../math_funcs.hpp"
7 #include "../parameter_input.hpp"
8 #include "../coordinates/coordinates.hpp"
9 #include "../thermodynamics/thermodynamics.hpp"
10 #include "../mesh/mesh.hpp"
11 #include "../utils/utils.hpp"
12 #include "../math/core.h"
13 #include "../debugger/debugger.hpp"
14 #include "radiation.hpp"
15 #include "radiation_utils.hpp" // setRadiationFlags
16 
17 Real const Radiation::hPlanck = 6.63E-34;
18 Real const Radiation::hPlanck_cgs = 6.63E-27;
19 Real const Radiation::cLight = 3.E8;
20 Real const Radiation::cLight_cgs = 3.E10;
21 Real const Radiation::stefanBoltzmann = 5.670374419E-8;
22 
24  pmy_block(pmb), pband(nullptr), rflags(0LL),
25  cooldown(0.), current(0.), planet(nullptr), stellarDistance_au_(1.)
26 {}
27 
29  pmy_block(pmb), pband(nullptr), rflags(0LL)
30 {
31  pmb->pdebug->Enter("Radiation");
32  std::stringstream &msg = pmb->pdebug->msg;
33 
34  // radiation flags
35  setRadiationFlags(&rflags, pin->GetOrAddString("radiation", "flags", ""));
36 
37  // distance to parent star
38  stellarDistance_au_ = pin->GetOrAddReal("radiation", "distance_au", 1.);
39  msg << "- stellar distance = " << stellarDistance_au_ << " au" << std::endl;
40 
41  // radiation bands
42  int bid = 1;
43  readRadiationBands(pin, bid);
44 
45  // incoming radiation direction (mu,phi) in degree
46  std::string str = pin->GetOrAddString("radiation", "indir", "(0.,0.)");
48 
49  // output radiance
50  radiance_units = pin->GetOrAddString("radiation", "radiance_units", "w/(m^2.cm^{-1}.sr)");
51  int nout = 0;
53  while (p != nullptr) {
54  nout += p->rayOutput.size();
55  p = p->next;
56  }
57  if (nout > 0) {
58  radiance.NewAthenaArray(nout, pmb->ncells3, pmb->ncells2);
59  // set band toa
60  p = pband;
61  nout = 0;
62  while (p != nullptr) {
63  p->btoa.InitWithShallowSlice(radiance,3,nout,p->rayOutput.size());
64  nout += p->rayOutput.size();
65  p = p->next;
66  }
67  }
68 
69  // time control
70  cooldown = pin->GetOrAddReal("radiation", "dt", 0.);
71  current = 0.;
72 
73  planet = new CelestrialBody(pin);
74  pmb->pdebug->Leave();
75 }
76 
78 {
79  if (pband != nullptr) {
80  while (pband->prev != nullptr) // should not be true
81  delete pband->prev;
82  while (pband->next != nullptr)
83  delete pband->next;
84  delete pband;
85  }
86  delete planet;
87 }
88 
90  std::stringstream msg;
92  int b = 0;
93  while (p != nullptr) {
94  if (b++ == n) break;
95  p = p->next;
96  }
97  return p;
98 }
99 
101  int n = 0;
102  RadiationBand* p = pband;
103  while (p != nullptr) {
104  p = p->next;
105  n++;
106  }
107  return n;
108 }
109 
111  int k, int j, int il, int iu)
112 {
113  pmy_block->pdebug->Call("Radiation::calculateRadiativeFluxes");
114  Coordinates *pcoord = pmy_block->pcoord;
115  Real dist = stellarDistance_au_;
116 
117  RadiationBand *p = pband;
118  if (pband == nullptr) return;
119 
120  while (p != nullptr) {
121  Direction ray;
122  if (p->bflags & RadiationFlags::Dynamic) {
123  planet->ParentZenithAngle(&ray.mu, &ray.phi, time, pcoord->x2v(j), pcoord->x3v(k));
124  dist = planet->ParentDistanceInAu(time);
125  } else {
126  ray = rayInput[0];
127  }
128 
129  // iu ~= ie + 1
130  p->setSpectralProperties(w, k, j, il - NGHOST, iu + NGHOST - 1);
131  p->calculateRadiativeFlux(ray, dist, k, j, il, iu);
132  p = p->next;
133  }
134  pmy_block->pdebug->Leave();
135 }
136 
138  int k, int j, int il, int iu)
139 {
140  pmy_block->pdebug->Call("Radiation::calculateRadiances");
141  Coordinates *pcoord = pmy_block->pcoord;
142  Real dist = stellarDistance_au_;
143 
144  RadiationBand *p = pband;
145  if (pband == nullptr) return;
146 
147  Direction ray;
149  planet->ParentZenithAngle(&ray.mu, &ray.phi, time, pcoord->x2v(j), pcoord->x3v(k));
150  dist = planet->ParentDistanceInAu(time);
151  }
152 
153  while (p != nullptr) {
154  // iu ~= ie + 1
155  p->setSpectralProperties(w, k, j, il - NGHOST, iu + NGHOST - 1);
156  p->calculateRadiance(ray, dist, k, j, il, iu);
157  p = p->next;
158  }
159  pmy_block->pdebug->Leave();
160 }
161 
163  int k, int j, int il, int iu)
164 {
165  RadiationBand *p = pband;
166  if (pband == nullptr) return;
167 
168  MeshBlock *pmb = pmy_block;
169  pmb->pdebug->Call("Radiation::addRadiativeFluxes");
170 
171  // x1-flux divergence
172  p = pband;
173  while (p != nullptr) {
174 #pragma omp simd
175  for (int i = il; i <= iu; ++i)
176  x1flux(IEN,k,j,i) += p->bflxup(k,j,i) - p->bflxdn(k,j,i);
177  p = p->next;
178  }
179  pmb->pdebug->Leave();
180 }
181 
183 {
184  char name[80];
185  RadiationBand *plast = pband;
186  while (plast != nullptr)
187  plast = plast->next;
188 
189  while (true) {
190  sprintf(name, "b%d", bid);
191  if (!pin->DoesParameterExist("radiation", name))
192  break;
193  RadiationBand* p = new RadiationBand(this, name, pin);
194  if (plast == nullptr) {
195  plast = p;
196  pband = p;
197  } else {
198  plast->next = p;
199  plast->next->prev = plast;
200  plast->next->next = nullptr;
201  plast = plast->next;
202  }
203  bid++;
204  }
205 
206  if (pin->DoesParameterExist("radiation", "bands_file")) {
207  ParameterInput* pin_next = new ParameterInput;
208  IOWrapper infile;
209  infile.Open(pin->GetString("radiation", "bands_file").c_str(), IOWrapper::FileMode::read);
210  pin_next->LoadFromFile(infile);
211  infile.Close();
212  InputBlock *pblock = pin->GetPtrToBlock("radiation");
213  InputLine *pline = pblock->pline;
214 
215  // remove the bands_file line
216  while (pline->pnext != nullptr) {
217  if (pline->pnext->param_name == "bands_file") {
218  InputLine *pnext = pline->pnext->pnext;
219  delete pline->pnext;
220  pline->pnext = pnext;
221  continue;
222  }
223  pline = pline->pnext;
224  }
225 
226  // get the first line of the current input in block radiation
227  pline = pin_next->GetPtrToBlock("radiation")->pline;
228 
229  // copy the current lines into the main input
230  while (pline != nullptr) {
231  pin->AddParameter(pblock, pline->param_name, pline->param_value, pline->param_comment);
232  pline = pline->pnext;
233  }
234  readRadiationBands(pin, bid);
235  delete pin_next;
236  }
237 }
238 
240  int num = 0;
241  RadiationBand *p = pband;
242  while (p != nullptr) {
243  num += p->rayOutput.size();
244  p = p->next;
245  }
246  return num;
247 }
248 
250 {
251  size_t size = 0;
252 
253  RadiationBand *p = pband;
254  while (p != nullptr) {
255  size += p->bflxup.GetSizeInBytes() + p->bflxdn.GetSizeInBytes();
256  p = p->next;
257  }
258 
259  return size;
260 }
261 
262 size_t Radiation::dumpRestartData(char *pdst)
263 {
264  RadiationBand *p = pband;
265  int offset = 0;
266  while (p != nullptr) {
267  std::memcpy(pdst + offset, p->bflxup.data(), p->bflxup.GetSizeInBytes());
268  offset += p->bflxup.GetSizeInBytes();
269  std::memcpy(pdst + offset, p->bflxdn.data(), p->bflxdn.GetSizeInBytes());
270  offset += p->bflxdn.GetSizeInBytes();
271  p = p->next;
272  }
273 
274  return getRestartDataSizeInBytes();
275 }
276 
277 size_t Radiation::loadRestartData(char *psrc)
278 {
279  RadiationBand *p = pband;
280  int offset = 0;
281  while (p != nullptr) {
282  std::memcpy(p->bflxup.data(), psrc + offset, p->bflxup.GetSizeInBytes());
283  offset += p->bflxup.GetSizeInBytes();
284  std::memcpy(p->bflxdn.data(), psrc + offset, p->bflxdn.GetSizeInBytes());
285  offset += p->bflxdn.GetSizeInBytes();
286  p = p->next;
287  }
288 
289  return getRestartDataSizeInBytes();
290 }
double Real
Definition: athena.hpp:29
@ IEN
Definition: athena.hpp:135
#define LL(iq, lc)
Definition: cdisort.h:293
AthenaArray< Real > x2v
Definition: coordinates.hpp:42
AthenaArray< Real > x3v
Definition: coordinates.hpp:42
InputLine * pline
Debugger * pdebug
Definition: mesh.hpp:140
int ncells3
Definition: mesh.hpp:97
Coordinates * pcoord
Definition: mesh.hpp:122
int ncells2
Definition: mesh.hpp:97
int DoesParameterExist(std::string block, std::string name)
void AddParameter(InputBlock *pib, std::string name, std::string value, std::string comment)
std::string GetString(std::string block, std::string name)
InputBlock * GetPtrToBlock(std::string name)
std::string GetOrAddString(std::string block, std::string name, std::string value)
void LoadFromFile(IOWrapper &input)
Real GetOrAddReal(std::string block, std::string name, Real value)
RadiationBand * prev
Definition: radiation.hpp:50
RadiationBand * next
Definition: radiation.hpp:50
static Real const cLight_cgs
Definition: radiation.hpp:97
void calculateRadiances(AthenaArray< Real > const &w, Real time, int k, int j, int il, int iu)
Definition: radiation.cpp:137
uint64_t rflags
Definition: radiation.hpp:103
MeshBlock * pmy_block
Definition: radiation.hpp:101
Real cooldown
Definition: radiation.hpp:104
static Real const hPlanck
Definition: radiation.hpp:94
size_t dumpRestartData(char *pdst)
Definition: radiation.cpp:262
Radiation(MeshBlock *pmb)
Definition: radiation.cpp:23
size_t loadRestartData(char *psrc)
Definition: radiation.cpp:277
Real stellarDistance_au_
Definition: radiation.hpp:134
void calculateRadiativeFluxes(AthenaArray< Real > const &w, Real time, int k, int j, int il, int iu)
Definition: radiation.cpp:110
static Real const hPlanck_cgs
Definition: radiation.hpp:95
int getNumBands()
Definition: radiation.cpp:100
static Real const stefanBoltzmann
Definition: radiation.hpp:98
std::string radiance_units
Definition: radiation.hpp:106
CelestrialBody * planet
Definition: radiation.hpp:105
RadiationBand * getBand(int n)
Definition: radiation.cpp:89
std::vector< Direction > rayInput
Definition: radiation.hpp:109
int getTotalNumberOutgoingRays()
Definition: radiation.cpp:239
size_t getRestartDataSizeInBytes()
Definition: radiation.cpp:249
void addRadiativeFluxes(AthenaArray< Real > &x1flux, int k, int j, int il, int iu)
Definition: radiation.cpp:162
void readRadiationBands(ParameterInput *pin, int &b)
Definition: radiation.cpp:182
RadiationBand * pband
Definition: radiation.hpp:102
Real current
Definition: radiation.hpp:104
static Real const cLight
Definition: radiation.hpp:96
AthenaArray< Real > radiance
Definition: radiation.hpp:110
const uint64_t Dynamic
Definition: radiation.hpp:34
void readRadiationDirections(std::vector< Direction > &ray, std::string str)
void setRadiationFlags(uint64_t *flags, std::string str)
Real phi
Definition: radiation.hpp:29
std::string param_value
std::string param_name
std::string param_comment
InputLine * pnext