cloudy trunk
Loading...
Searching...
No Matches
conv_base.cpp
Go to the documentation of this file.
1/* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
2 * others. For conditions of distribution and use see copyright notice in license.txt */
3/*ConvBase main routine to drive ionization solution for all species, find total opacity
4 * called by ConvIoniz */
5/*lgConverg check whether ionization of element nelem has converged */
6#include "cddefines.h"
7#include "dynamics.h"
8#include "trace.h"
9#include "elementnames.h"
10#include "save.h"
11#include "phycon.h"
12#include "secondaries.h"
13#include "stopcalc.h"
14#include "grainvar.h"
15#include "highen.h"
16#include "dense.h"
17#include "hmi.h"
18#include "rfield.h"
19#include "pressure.h"
20#include "taulines.h"
21#include "rt.h"
22#include "grains.h"
23#include "atmdat.h"
24#include "ionbal.h"
25#include "opacity.h"
26#include "cooling.h"
27#include "thermal.h"
28#include "mole.h"
29#include "iso.h"
30#include "conv.h"
31#include "h2.h"
32#include "deuterium.h"
33
34STATIC bool lgNetEdenSrcSmall( void );
35
36void UpdateUTAs( void );
37
38/*lgIonizConverg check whether ionization of element nelem has converged, called by ionize,
39 * returns true if element is converged, false if not */
40namespace
41{
42 class IonizConverg
43 {
44 /* this is fractions [ion stage][nelem], ion stage = 0 for atom, nelem=0 for H*/
45 double OldFracs[LIMELM+1][LIMELM+1];
46 public:
47 IonizConverg()
48 {
49 for (long nelem = ipHYDROGEN; nelem < LIMELM; ++nelem)
50 {
51 for (long ion = 0; ion <= nelem+1; ++ion)
52 {
53 OldFracs[nelem][ion] = dense.xIonDense[nelem][ion];
54 }
55 }
56 }
57 void operator()(
58 /* atomic number on the C scale, 0 for H, 25 for Fe */
59 long loop_ion ,
60 /* this is allowed error as a fractional change. Most are 0.15 */
61 double delta ,
62 /* option to print abundances */
63 bool lgPrint )
64 {
65 bool lgConverg_v = false;
66 double Abund,
67 bigchange ,
68 change ,
69 one;
70
71 DEBUG_ENTRY( "IonizConverg::operator()" );
72
73 for (long nelem =ipHYDROGEN; nelem<LIMELM; ++nelem )
74 {
75 if( !dense.lgElmtOn[nelem] )
76 continue;
77
78 double abundold=0. , abundnew=0.;
79 long ionchg=-1;
80
81 /* arguments are atomic number, ionization stage
82 * OldFracs[nelem][0] is abundance of nelem (cm^-3) */
83
84 /* this function returns true if ionization of element
85 * with atomic number nelem has not changed by more than delta,*/
86
87 /* check whether abundances exist yet, only do this for after first zone */
88 /*if( nzone > 0 )*/
89 /* >>chng 03 sep 02, check on changed ionization after first call through,
90 * to insure converged constant eden / temperature models */
91 if( conv.nPres2Ioniz )
92 {
93 /* >>chng 04 aug 31, this had been static, caused last hits on large changes
94 * in ionization */
95 bigchange = 0.;
96 change = 0.;
97 Abund = dense.gas_phase[nelem];
98
99 /* loop over all ionization stages, loop over all ions, not just active ones,
100 * since this also sets old values, and so will zero out non-existant ones */
101 for( long ion=0; ion <= (nelem+1); ++ion )
102 {
103 /*lint -e727 OlsdFracs not initialized */
104 if( OldFracs[nelem][ion]/Abund > 1e-4 &&
105 dense.xIonDense[nelem][ion]/Abund > 1e-4 )
106 /*lint +e727 OlsdFracs not initialized */
107 {
108 /* check change in old to new value */
109 one = fabs(dense.xIonDense[nelem][ion]-OldFracs[nelem][ion])/
110 OldFracs[nelem][ion];
111 change = MAX2(change, one );
112 /* remember abundances for largest change */
113 if( change>bigchange )
114 {
115 bigchange = change;
116 abundold = OldFracs[nelem][ion]/Abund;
117 abundnew = dense.xIonDense[nelem][ion]/Abund;
118 ionchg = ion;
119 }
120 }
121 /* now set new value */
122 OldFracs[nelem][ion] = dense.xIonDense[nelem][ion];
123 }
124
125 if( change >= delta )
126 {
127 char chConvIoniz[INPUT_LINE_LENGTH];
128 sprintf( chConvIoniz , "%2s ion" , elementnames.chElementSym[nelem] );
129 conv.setConvIonizFail(chConvIoniz,abundold,abundnew);
130 ASSERT( abundold>0. && abundnew>0. );
131 }
132 else
133 lgConverg_v = true;
134 }
135 else
136 {
137 for( long ion=0; ion <= (nelem+1); ++ion )
138 {
139 OldFracs[nelem][ion] = dense.xIonDense[nelem][ion];
140 }
141
142 }
143
144 /* option to print abundances */
145 if( lgPrint )
146 {
147 fprintf(ioQQQ," nz %ld loop %ld element %li converged? %c worst %ld change %g\n",
148 nzone, loop_ion, nelem, TorF(lgConverg_v),ionchg,bigchange);
149 for( long ion=0; ion<(nelem+1); ++ion )
150 {
151 fprintf(ioQQQ,"\t%.5e", dense.xIonDense[nelem][ion]/dense.gas_phase[nelem]);
152 }
153 fprintf(ioQQQ,"\n");
154 }
155 }
156 }
157 };
158}
159
160/*ConvBase main routine to drive ionization solution for all species, find total opacity
161 * called by ConvIoniz
162 * return 0 if ok, 1 if abort */
164 /* this is zero the first time ConvBase is called by convIoniz,
165 * counts number of call thereafter */
166 long loopi )
167{
168 double HeatOld,
169 EdenTrue_old,
170 EdenFromMolecOld,
171 EdenFromGrainsOld,
172 HeatingOld ,
173 CoolingOld;
174 static double SecondOld;
175 static long int nzoneOTS=-1;
176 const int LOOP_ION_LIMIT = 10;
177 long int loop_ion;
178 static double SumOTS=0. , OldSumOTS[2]={0.,0.};
179 double save_iso_grnd[NISO][LIMELM];
180 valarray<realnum> mole_save(mole_global.num_calc);
181 IonizConverg lgIonizConverg;
182
183 DEBUG_ENTRY( "ConvBase()" );
184
185 /* this is set to phycon.te in tfidle, is used to insure that all temp
186 * vars are properly updated when conv_ionizeopacitydo is called
187 * NB must be same type as phycon.te */
188 ASSERT( fp_equal( phycon.te, thermal.te_update ) );
189
190 /* this allows zone number to be printed with slight increment as zone converged
191 * conv.nPres2Ioniz is incremented at the bottom of this routine */
192 fnzone = (double)nzone + (double)conv.nPres2Ioniz/100.;
193
194 /* reevaluate pressure */
195 /* this sets values of pressure.PresTotlCurr, also calls tfidle,
196 * and sets the total energy content of gas, which may be important for acvection */
198
199 /* >>chng 04 sep 15, move EdenTrue_old up here, and will redo at bottom
200 * to find change
201 * find and save current true electron density - if this changes by more than the
202 * tolerance then ionization solution is not converged */
203 /* >>chng 04 jul 27, move eden_sum here from after this loop, so that change in EdenTrue
204 * can be monitored */
205 /* update EdenTrue, eden itself is actually changed in ConvEdenIoniz */
206 /* must not call eden_sum on very first time since for classic PDR total
207 * ionization may still be zero on first call */
208 if( conv.nTotalIoniz )
209 {
210 if( eden_sum() )
211 {
212 /* non-zero return indicates abort condition */
213 ++conv.nTotalIoniz;
214 return 1;
215 }
216 }
217
218 /* the following two were evaluated in eden_sum
219 * will confirm that these are converged */
220 EdenTrue_old = dense.EdenTrue;
221 EdenFromMolecOld = mole.elec;
222 EdenFromGrainsOld = gv.TotalEden;
223 HeatingOld = thermal.htot;
224 CoolingOld = thermal.ctot;
225
226 /* remember current ground state population - will check if converged */
227 for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
228 {
229 for( long nelem=ipISO; nelem<LIMELM;++nelem )
230 {
231 if( dense.lgElmtOn[nelem] )
232 {
233 /* save the ground state population */
234 save_iso_grnd[ipISO][nelem] = iso_sp[ipISO][nelem].st[0].Pop();
235 }
236 }
237 }
238
239 for( long i=0; i < mole_global.num_calc; i++ )
240 {
241 mole_save[i] = (realnum) mole.species[i].den;
242 }
243
244 if( loopi==0 )
245 {
246 /* these will be used to look for oscillating ots rates */
247 OldSumOTS[0] = 0.;
248 OldSumOTS[1] = 0.;
249 conv.lgOscilOTS = false;
250 }
251
252 if( trace.lgTrace )
253 {
254 fprintf( ioQQQ,
255 " ConvBase called. %.2f Te:%.3e HI:%.3e HII:%.3e H2:%.3e Ne:%.3e htot:%.3e CSUP:%.2e Conv?%c\n",
256 fnzone,
257 phycon.te,
258 dense.xIonDense[ipHYDROGEN][0],
259 dense.xIonDense[ipHYDROGEN][1],
260 hmi.H2_total,
261 dense.eden,
262 thermal.htot,
263 secondaries.csupra[ipHYDROGEN][0] ,
264 TorF(conv.lgConvIoniz()) );
265 }
266 /* want this flag to be true when we exit, various problems will set falst */
267 conv.resetConvIoniz();
268
269 /* this routine is in heatsum.c, and zeros out the heating array */
270 HeatZero();
271
272 /* if this is very first call, say not converged, since no meaningful way to
273 * check on changes in quantities. this counter is false only on very first
274 * call, when equal to zero */
275 if( !conv.nTotalIoniz )
276 {
277 conv.setConvIonizFail( "first call", 0.0, 0.0 );
278 }
279
280 /* this will be flag to check whether any ionization stages
281 * were trimmed down */
282 conv.lgIonStageTrimed = false;
283
284 /* must redo photoionization rates for all species on second and later tries */
285 /* always reevaluate the rates when . . . */
286 /* this flag says not to redo opac and photo rates, and following test will set
287 * true if one of several tests not done*/
288 opac.lgRedoStatic = false;
289 if(
290 /* opac.lgOpacStatic (usually true), is set false with no static opacity command */
291 !opac.lgOpacStatic ||
292 /* we are in search mode */
293 conv.lgSearch ||
294 /* this is the first call to this zone */
295 conv.nPres2Ioniz == 0 )
296 {
297 /* we need to redo ALL photoionization rates */
298 opac.lgRedoStatic = true;
299 }
300
301 /* calculate radiative and dielectronic recombination rate coefficients */
303
304 /* this adjusts the lowest and highest stages of ionization we will consider,
305 * only safe to call when lgRedoStatic is true since this could lower the
306 * lowest stage of ionization, which needs all its photo rates */
307
308 /* conv.nTotalIoniz is only 0 (false) on the very first call to here,
309 * when the ionization distribution is not yet done */
310 if( conv.nTotalIoniz )
311 {
312 bool lgIonizTrimCalled = false;
313 static long int nZoneCalled = 0;
314
315 fixit(); // nZoneCalled should be reinitialized for each grid point?
316
317 /* ionization trimming only used for He and heavier, not H */
318 /* only do this one time per zone since up and down cycle can occur */
319 /* >>chng 05 jan 15, increasing temperature above default first conditions, also
320 * no trim during search - this fixed major logical error when sim is
321 * totally mechanically heated to coronal temperatures -
322 * problem discovered by Ronnie Hoogerwerf */
323 /* do not keep changing the trim after the first call within
324 * this zone once we are deep into layer - doing so introduced very
325 * small level of noise as some stages
326 * went up and down - O+2 - O+3 was good example, when small H+3 after He+ i-front
327 * limit to one increase per element per zone */
328 if( conv.nTotalIoniz>2 &&
329 /* only call one time per zone except during search phase,
330 * when only call after 20 times only if temperature has changed */
331 ( conv.lgSearch || nZoneCalled!=nzone) )
332 {
333 lgIonizTrimCalled = true;
334 for( long nelem=ipHELIUM; nelem<LIMELM; ++nelem )
335 {
336 if( dense.lgElmtOn[nelem] )
337 {
338 /* ion_trim will set conv.lgIonStageTrimed true is any ion has its
339 * lowest stage of ionization dropped or trimmed */
340 ion_trim(nelem);
341 }
342 }
343 nZoneCalled = nzone;
344 }
345
346 /* following block only set of asserts */
347# if !defined(NDEBUG)
348 /* check that proper abundances are either positive or zero */
349 for( long nelem=ipHELIUM; nelem<LIMELM; ++nelem)
350 {
351 if( dense.lgElmtOn[nelem] )
352 {
353 for( long ion=0; ion<dense.IonLow[nelem]; ++ion )
354 {
355 ASSERT( dense.xIonDense[nelem][ion] == 0. );
356 }
357 /*if( nelem==5 ) fprintf(ioQQQ,"carbbb\t%li\n", dense.IonHigh[nelem]);*/
358 for( long ion=dense.IonLow[nelem]; ion<=dense.IonHigh[nelem]; ++ion )
359 {
360 /* >>chng 02 feb 06, had been > o., chng to > SMALLFLOAT to
361 * trip over VERY small floats that failed on alphas, but not 386
362 *
363 * in case where lower ionization stage was just lowered or
364 * trimmed down the abundance
365 * was set to SMALLFLOAT so test must be < SMALLFLOAT */
366 /* >>chng 02 feb 19, add check for search phase. During this search
367 * models with extreme ionization (all neutral or all ionized) can
368 * have extreme but non-zero abundances far from the ionization peak for
369 * element with lots of electrons. These will go away once the model
370 * becomes stable */
371 /* >>chng 03 dec 01, add check on whether ion trim was called
372 * conserve.in threw assert when iontrim not called and abund grew small */
373 ASSERT( conv.lgSearch || !lgIonizTrimCalled ||
374 /* this can happen if all C is in the form of CO */
375 (ion==0 && dense.IonHigh[nelem]==0 ) ||
376 dense.xIonDense[nelem][ion] >= SMALLFLOAT ||
377 dense.xIonDense[nelem][ion]/dense.gas_phase[nelem] >= SMALLFLOAT );
378 }
379 for( long ion=dense.IonHigh[nelem]+1; ion<nelem+1; ++ion )
380 {
381 ASSERT( ion >= 0 );
382 ASSERT( dense.xIonDense[nelem][ion] == 0. );
383 }
384 }
385 }
386# endif
387 }
388
389 /* now check if anything trimmed down */
390 if( conv.lgIonStageTrimed )
391 {
392 /* something was trimmed down, so say that ionization not yet stable */
393 /* say that ionization has not converged, secondaries changed by too much */
394 conv.setConvIonizFail( "IonTrimmed", 0.0, 0.0 );
395 }
396
397 /* reevaluate advective terms if turned on */
398 if( dynamics.lgAdvection )
399 DynaIonize();
400
401 /* evaluate Compton heating, bound E Compton, cosmic rays */
402 highen();
403
404 // Depends on dense.eden, phycon.te and trimming level of H, He
406
408
409 for( diatom_iter diatom = diatoms.begin(); diatom != diatoms.end(); ++diatom )
410 {
411 (*diatom)->CalcPhotoionizationRate();
412 }
413
414 long ionLowCache[LIMELM], ionHighCache[LIMELM];
415 for( long nelem=ipHYDROGEN; nelem<LIMELM; nelem++ )
416 {
417 if( dense.lgSetIoniz[nelem] )
418 {
419 dense.IonLow[nelem] = 0;
420 dense.IonHigh[nelem] = nelem + 1;
421 while( dense.SetIoniz[nelem][dense.IonLow[nelem]] == 0. )
422 ++dense.IonLow[nelem];
423 while( dense.SetIoniz[nelem][dense.IonHigh[nelem]] == 0. )
424 --dense.IonHigh[nelem];
425 }
426 ionLowCache[nelem] = dense.IonLow[nelem];
427 ionHighCache[nelem] = dense.IonHigh[nelem];
428 }
429
430 /* >>chng 04 feb 15, add loop over ionization until converged. Non-convergence
431 * in this case means ionization ladder pop changed, probably because of way
432 * that Auger ejection is treated - loop exits when conditions tested at end are met */
433 conv.incrementCounter(CONV_BASE_CALLS);
434 loop_ion = 0;
435 do
436 {
437 conv.incrementCounter(CONV_BASE_LOOPS);
439
440 /* set the convergence flag to true,
441 * if anything changes in ConvBase, it will be set false */
442 if( loop_ion )
443 conv.resetConvIoniz();
444
445 /* charge transfer evaluation needs to be here so that same rate
446 * coefficient used for H ion and other recombination */
447 /* fill in master charge transfer array, and zero out some arrays that track effects */
448
450
451 /* find grain temperature, charge, and gas heating rate */
452 /* >>chng 04 apr 15, moved here from after call to HeLike(), PvH */
453 GrainDrive();
454
455 /* evaluate Compton heating, bound E Compton, cosmic rays */
456 highen();
457
458 /* find corrections for three-body rec - collisional ionization */
459 atmdat_3body();
460
461 /* update UTA inner-shell ionization rates */
462 UpdateUTAs();
463
465
466 long ion_loop = 0;
467 const int nconv = 3;
468 double xIonDense0[nconv][LIMELM][LIMELM+1];
469 bool lgShortCircuit = false;
470 for( ion_loop=0; ion_loop<nconv && !lgShortCircuit; ++ion_loop)
471 {
472 conv.incrementCounter(CONV_BASE_ACCELS);
473 for (long nelem = ipHYDROGEN; nelem < LIMELM; ++nelem)
474 {
475 for (long ion = 0; ion <= nelem+1; ++ion)
476 {
477 xIonDense0[ion_loop][nelem][ion] = dense.xIonDense[nelem][ion];
478 }
479 }
480
481 double netion = 0.0;
482 for( long nelem=ipHYDROGEN; nelem<LIMELM; nelem++ )
483 {
484 if (!dense.lgElmtOn[nelem])
485 continue;
486
487 ASSERT(dense.IonLow[nelem] >= ionLowCache[nelem]);
488 ASSERT(dense.IonHigh[nelem] <= ionHighCache[nelem]);
489 ion_wrapper( nelem );
490 if (nelem == ipHYDROGEN || nelem == ipOXYGEN)
491 {
492 double hion =
493 atmdat.CharExcRecTo[ipHYDROGEN][ipOXYGEN][0]*
494 iso_sp[ipH_LIKE][ipHYDROGEN].st[0].Pop()*
495 dense.xIonDense[ipOXYGEN][1]-
496 atmdat.CharExcIonOf[ipHYDROGEN][ipOXYGEN][0]*
497 dense.xIonDense[ipHYDROGEN][1]*
498 dense.xIonDense[ipOXYGEN][0];
499 if (nelem == ipHYDROGEN)
500 netion += hion;
501 else
502 netion -= hion;
503 }
504
505 if ( conv.lgUpdateCouplings && conv.nTotalIoniz )
507 }
508 if (dense.lgElmtOn[ipHYDROGEN] && dense.lgElmtOn[ipOXYGEN] && ion_loop == 1)
509 {
510 // Comparison rate is whether error is either a small fraction of the minimum ionization rate,
511 // or a tighter check on whether the fluxes are balanced as well as is feasible numerically
512 double ion_cmp = MAX2(0.01*MIN2(ionbal.RateIonizTot(ipHYDROGEN,0), ionbal.RateIonizTot(ipOXYGEN,0)),
513 1e-4*atmdat.CharExcRecTo[ipHYDROGEN][ipOXYGEN][0]*
514 iso_sp[ipH_LIKE][ipHYDROGEN].st[0].Pop()*dense.xIonDense[ipOXYGEN][1]);
515
516 if (fabs(netion) > ion_cmp &&
517 dense.xIonDense[ipOXYGEN][1]*dense.xIonDense[ipHYDROGEN][1] >
518 1e-8*iso_sp[ipH_LIKE][ipHYDROGEN].st[0].Pop()*dense.xIonDense[ipOXYGEN][0] )
519 {
520 conv.setConvIonizFail( "OH CX inconsistency" , ion_cmp, netion);
521 }
522 }
523
524 // If not going to end anyhow, check whether changes were small
525 bool lgCanShortCircuit = (ion_loop+1 < nconv);
526 for( long nelem=ipHYDROGEN; nelem<LIMELM && lgCanShortCircuit; ++nelem )
527 {
528 if (!dense.lgElmtOn[nelem])
529 continue;
530 for (long ion = dense.IonLow[nelem];
531 ion <= dense.IonHigh[nelem]; ++ion)
532 {
533 double x0 = xIonDense0[ion_loop][nelem][ion];
534 double x1 = dense.xIonDense[nelem][ion];
535 if (fabs(x0-x1) > 1e-6*(x0+x1))
536 {
537 lgCanShortCircuit = false;
538 break;
539 }
540 }
541 }
542 lgShortCircuit = lgCanShortCircuit;
543 }
544
545 if (!lgShortCircuit)
546 {
547 // Apply convergence acceleration
548 for (long nelem = ipHYDROGEN; nelem < LIMELM; ++nelem)
549 {
550 if (!dense.lgElmtOn[nelem])
551 continue;
552 double tot0 = 0., tot1 = 0.;
553 double xIonNew[LIMELM+1];
554 for (long ion = 0; ion <= nelem+1; ++ion)
555 {
556 double x0 = xIonDense0[nconv-2][nelem][ion];
557 double x1 = xIonDense0[nconv-1][nelem][ion];
558 double x2 = dense.xIonDense[nelem][ion];
559 xIonNew[ion] = x2;
560 tot0 += x2;
561 // Richardson extrapolation formula to accelerate convergence
562 // Assumes convergence to x^ is geometric, i.e.
563 // x_i = x^ + a d^i
564
565 double extstep = 0.,predict=x2,
566 step0 = x1-x0, step1 = x2-x1, abs1 = fabs(step1);
567 // Protect against roundoff/noise in inner solver
568 if ( abs1 > 1000.0*((double)DBL_EPSILON)*x2 )
569 {
570 double denom = fabs(step1-step0);
571 double sgn = (step1*step0 > 0)? 1.0 : -1.0;
572 // Greatest acceleration allowed is MAXACC*latest step length
573 // Can we do better than static tuning this parameter?
574 const double MAXACC=100.0;
575 double extfac = 1.0/(denom/abs1 + 1.0/MAXACC);
576 extstep = sgn*extfac*step1;
577 //extstep = sgn*MIN2(extfac*step1,0.01*x2);
578 predict = x2+extstep;
579 if (predict > 0.0)
580 xIonNew[ion] = predict;
581 }
582 if ( 0 )
583 //if ( nelem == ipHYDROGEN || (nelem == ipIRON && ion <=2 ) )
584 if ( (nelem == ipNICKEL && ion <=2 ) )
585 fprintf(ioQQQ,"Extrap %3ld %3ld %13.6g %13.6g %13.6g %13.6g %13.6g %13.6g\n",
586 nelem,ion,
587 x0,x0-xIonDense0[nconv-3][nelem][ion],x1-x0,x2-x1,extstep,predict);
588 tot1 += xIonNew[ion];
589 }
590 if ( tot1 > SMALLFLOAT )
591 {
592 double scal = tot0/tot1;
593 for (long ion = 0; ion <= nelem+1; ++ion)
594 {
595 dense.xIonDense[nelem][ion] = scal*xIonNew[ion];
596 }
597 for ( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
598 {
599 double renorm;
600 iso_renorm(nelem, ipISO, renorm);
601 }
602 }
603 }
604
605 bool lgPostExtrapSolve = true;
606 if (lgPostExtrapSolve)
607 {
608 for (long nelem = ipHYDROGEN; nelem < LIMELM; ++nelem)
609 {
610 if (!dense.lgElmtOn[nelem])
611 continue;
612 ion_wrapper( nelem );
613 }
614 }
615 }
616 SetDeuteriumIonization( dense.xIonDense[ipHYDROGEN][0], dense.xIonDense[ipHYDROGEN][1] );
617
618 /*>>chng 04 may 09, add option to abort here, inspired by H2 pop failures
619 * which cascaded downstream if we did not abort */
620 /* return if one of above solvers has declared abort condition */
621 if( lgAbort )
622 {
623 ++conv.nTotalIoniz;
624 return 1;
625 }
626
628
629 /* drive chemistry network*/
630 mole_drive();
631
633
634 /* all elements have now had ionization reevaluated - in some cases we may have upset
635 * the ionization that was forced with an "element ionization" command - here we will
636 * re-impose that set ionization */
637 /* >>chng 04 oct 13, add this logic */
638 for( long nelem=ipHYDROGEN; nelem < LIMELM; nelem++ )
639 {
640 if( dense.lgSetIoniz[nelem] )
641 {
642 dense.IonLow[nelem] = 0;
643 dense.IonHigh[nelem] = nelem + 1;
644 while( dense.SetIoniz[nelem][dense.IonLow[nelem]] == 0. )
645 ++dense.IonLow[nelem];
646 while( dense.SetIoniz[nelem][dense.IonHigh[nelem]] == 0. )
647 --dense.IonHigh[nelem];
648 }
649 }
650
651 /* redo ct ion rate for reasons now unclear */
653
654 /* evaluate molecular hydrogen level populations */
655 for( diatom_iter diatom = diatoms.begin(); diatom != diatoms.end(); ++diatom )
656 {
657 bool lgPopsConverged = true;
658 double old_val, new_val;
659 (*diatom)->H2_LevelPops( lgPopsConverged, old_val, new_val );
660 if( !lgPopsConverged )
661 {
662 conv.setConvIonizFail( "H2 pops", old_val, new_val);
663 }
664 }
665
667
668 /* lgIonizConverg is a function to check whether ionization has converged
669 * check whether ionization changed by more than relative error
670 * given by second number */
671 /* >>chng 04 feb 14, loop over all elements rather than just a few */
672 lgIonizConverg(loop_ion, conv.IonizErrorAllowed , false );
673
674 if( deut.lgElmtOn )
675 {
676 static double OldDeut[2] = {0., 0.};
677 for( long ion=0; ion<2; ++ion )
678 {
679 if( fabs(deut.xIonDense[ion] - OldDeut[ion] ) > 0.2*conv.IonizErrorAllowed*fabs(OldDeut[ion]) )
680 {
681 conv.setConvIonizFail( "D ion" , OldDeut[ion], deut.xIonDense[ion]);
682 }
683 OldDeut[ion] = deut.xIonDense[ion];
684 }
685 }
686
687 if( fabs(EdenTrue_old - dense.EdenTrue) > conv.EdenErrorAllowed/2.*fabs(dense.EdenTrue) )
688 {
689 conv.setConvIonizFail( "EdTr cng" , EdenTrue_old, dense.EdenTrue);
690 }
691
692 for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
693 {
694 for( long nelem=ipISO; nelem<LIMELM; ++nelem )
695 {
696 lgStatesConserved( nelem, nelem-ipISO, iso_sp[ipISO][nelem].st, iso_sp[ipISO][nelem].numLevels_local, 1e-3f, loop_ion );
697 }
698 }
699
700 if( trace.nTrConvg>=4 )
701 {
702 /* trace ionization */
703 fprintf( ioQQQ,
704 " ConvBase4 ionization driver loop_ion %li converged? %c reason not converged %s\n" ,
705 loop_ion ,
706 TorF( conv.lgConvIoniz()) ,
707 conv.chConvIoniz() );
708 }
709
710 ++loop_ion;
711 }
712 /* loop is not converged, less than loop limit, and we are reevaluating */
713 while( !conv.lgConvIoniz() && loop_ion < LOOP_ION_LIMIT && rfield.lgIonizReevaluate);
714
715 if( conv.lgConvIoniz() )
717
718 /* >>chng 05 oct 29, move CT heating here from heat_sum since this sometimes contributes
719 * cooling rather than heat and so needs to be sorted out before either heating or cooling
720 * are derived first find net heating - */
721 thermal.char_tran_heat = ChargTranSumHeat();
722 /* net is cooling if negative */
723 thermal.char_tran_cool = MAX2(0. , -thermal.char_tran_heat );
724 thermal.char_tran_heat = MAX2(0. , thermal.char_tran_heat );
725
726 /* get total cooling, thermal.ctot = does not occur since passes as pointer. This can add heat.
727 * it calls coolSum at end to sum up the total cooling */
728 CoolEvaluate(&thermal.ctot );
729
730 /* get total heating rate - first save old quantities to check how much it changes */
731 HeatOld = thermal.htot;
732
733 /* HeatSum will update ElecFrac,
734 * secondary electron ionization and excitation efficiencies,
735 * and sum up current secondary rates - remember old values before we enter */
736 SecondOld = secondaries.csupra[ipHYDROGEN][0];
737
738 /* update the total heating - it was all set to zero in HeatZero at top of this routine,
739 * occurs before secondaries bit below, since this updates electron fracs */
740 HeatSum();
741
742 /* test whether we can just set the rate to the new one, or whether we should
743 * take average and do it again. secondaries.sec2total was set in hydrogenic, and
744 * is ratio of secondary to total hydrogen destruction rates */
745 /* >>chng 02 nov 20, add test on size of csupra - primal had very close to underflow */
746 if( (secondaries.csupra[ipHYDROGEN][0]>SMALLFLOAT) && (secondaries.sec2total>0.001) &&
747 fabs( 1. - SecondOld/SDIV(secondaries.csupra[ipHYDROGEN][0]) ) > 0.1 &&
748 SecondOld > 0. && secondaries.csupra[ipHYDROGEN][0] > 0.)
749 {
750 /* say that ionization has not converged, secondaries changed by too much */
751 conv.setConvIonizFail( "SecIonRate", SecondOld,
752 secondaries.csupra[ipHYDROGEN][0] );
753 }
754
755# if 0
756 static realnum hminus_old=0.;
757 /* >>chng 04 apr 15, add this convergence test */
758 if( conv.nTotalIoniz )
759 {
760 realnum hminus_den = findspecieslocal("H-")->den;
761 if( fabs( hminus_old-hminus_den ) > fabs( hminus_den * conv.EdenErrorAllowed ) )
762 {
763 conv.setConvIonizFail( "Big H- chn", hminus_old, hminus_den );
764 }
765 hminus_old = hminus_den;
766 }
767# endif
768
769 if( HeatOld > 0. && !thermal.lgTemperatureConstant )
770 {
771 /* check if heating has converged - tolerance is final match */
772 if( fabs(1.-thermal.htot/HeatOld) > conv.HeatCoolRelErrorAllowed*.5 )
773 {
774 conv.setConvIonizFail( "Big d Heat", HeatOld, thermal.htot);
775 }
776 }
777
778 /* abort flag may have already been set - if so bail */
779 if( lgAbort )
780 {
781
782 return 1;
783 }
784
785 /* evaluate current opacities
786 * rfield.lgOpacityReevaluate normally true,
787 * set false with no opacity reevaluate command, an option to only
788 * evaluate opacity one time per zone */
789 if( !conv.nPres2Ioniz || rfield.lgOpacityReevaluate || conv.lgIonStageTrimed )
791
792 /* >>chng 02 jun 11, call even first time that this routine is called -
793 * this seemed to help convergence */
794
795 /* do OTS rates for all lines and all continua since
796 * we now have ionization balance of all species. Note that this is not
797 * entirely self-consistent, since destruction probabilities here are not the same as
798 * the ones used in the model atoms. Problems?? if near convergence
799 * then should be nearly identical */
800 if( !conv.nPres2Ioniz || rfield.lgOpacityReevaluate || rfield.lgIonizReevaluate ||
801 conv.lgIonStageTrimed || conv.lgSearch || nzone!=nzoneOTS )
802 {
803 RT_OTS();
804 nzoneOTS = nzone;
805
806 /* remember old ots rates */
807 OldSumOTS[0] = OldSumOTS[1];
808 OldSumOTS[1] = SumOTS;
809 /*fprintf(ioQQQ," calling RT_OTS zone %.2f SumOTS is %.2e\n",fnzone,SumOTS);*/
810
811 /* now update several components of the continuum, this only happens after
812 * we have gone through entire solution for this zone at least one time.
813 * there can be wild ots oscillation on first call */
814 /* the rel change of 0.2 was chosen by running hizqso - smaller increased
815 * itrzn but larger did not further decrease it. */
816 RT_OTS_Update(&SumOTS);
817 /*fprintf(ioQQQ,"RT_OTS_Updateee\t%.3f\t%.2e\t%.2e\n", fnzone,SumOTS , OldSumOTS[1] );*/
818 }
819 else
820 SumOTS = 0.;
821
822 /* now check whether the ots rates changed */
823 if( SumOTS> 0. )
824 {
825 /* the ots rate must be converged to the error in the electron density */
826 /* how fine should this be converged?? originally had above, 10%, but take
827 * smaller ratio?? */
828 if( fabs(1.-OldSumOTS[1]/SumOTS) > conv.EdenErrorAllowed )
829 {
830 /* this branch, ionization not converged due to large change in ots rates.
831 * check whether ots rates are oscillating, if loopi > 1 so we have enough info*/
832 if( loopi > 1 )
833 {
834 /* here we have three values, are they changing sign? */
835 if( (OldSumOTS[0]-OldSumOTS[1]) * ( OldSumOTS[1] - SumOTS ) < 0. )
836 {
837 /* ots rates are oscillating */
838 conv.lgOscilOTS = true;
839 }
840 }
841
842 conv.setConvIonizFail( "OTSRatChng" , OldSumOTS[1], SumOTS);
843 }
844
845 /* produce info on the ots fields if either "trace ots" or
846 * "trace convergence xxx ots " was entered */
847 if( ( trace.lgTrace && trace.lgOTSBug ) ||
848 ( trace.nTrConvg && trace.lgOTSBug ) )
849 {
850 RT_OTS_PrtRate(SumOTS*0.05 , 'b' );
851 }
852 /*fprintf(ioQQQ,"DEBUG opac\t%.2f\t%.3e\t%.3e\n",fnzone,
853 dense.xIonDense[ipNICKEL][0] ,
854 dense.xIonDense[ipZINC][0] );*/
855 {
856 /* DEBUG OTS rates - turn this on to debug line, continuum or both rates */
857 enum {DEBUG_LOC=false};
858 if( DEBUG_LOC && (nzone>110)/**/ )
859 {
860# if 0
861# include "lines_service.h"
862 DumpLine( &iso_sp[ipH_LIKE][ipHYDROGEN].trans(2,0) );
863# endif
864 /* last par 'l' for line, 'c' for continua, 'b' for both,
865 * the numbers printed are:
866 * cell i, [i], so 1 less than ipoint
867 * anu[i],
868 * otslin[i],
869 * opacity_abs[i],
870 * otslin[i]*opacity_abs[i],
871 * rfield.chLineLabel[i] ,
872 * rfield.line_count[i] */
873 }
874 }
875 }
876 {
877 /* DEBUG OTS rates - turn this on to debug line, continuum or both rates */
878 enum {DEBUG_LOC=false};
879 if( DEBUG_LOC && (nzone>200)/**/ )
880 {
881 fprintf(ioQQQ,"debug otsss\t%li\t%.3e\t%.3e\t%.3e\n",
882 nzone,
883 iso_sp[0][1].trans(15,3).Emis().ots(),
884 TauLines[26].Emis().ots(),
885 opac.opacity_abs[2069]);
886 }
887 }
888
889 /* option to print OTS continuum with TRACE OTS */
890 if( trace.lgTrace && trace.lgOTSBug )
891 {
892 /* find ots rates here, so we only print fraction of it,
893 * SumOTS is both line and continuum contributing to ots, and is multiplied by opacity */
894 /* number is weakest rate to print */
895 RT_OTS_PrtRate(SumOTS*0.05 , 'b' );
896 }
897
898 // all RT routines called
899 RT_line_all( );
900
901 /* >>chng 01 mar 16, evaluate pressure here since changing and other values needed */
902 /* reevaluate pressure */
903 /* this sets values of pressure.PresTotlCurr, also calls tfidle */
905
907
908 /* update some counters that keep track of how many times this routine
909 * has been called */
910
911 if( trace.lgTrace )
912 {
913 fprintf( ioQQQ,
914 " ConvBase return. fnzone %.2f nPres2Ioniz %li Te:%.3e HI:%.3e HII:%.3e H2:%.3e Ne:%.3e htot:%.3e CSUP:%.2e Conv?%c reason:%s\n",
915 fnzone,
916 conv.nPres2Ioniz,
917 phycon.te,
918 dense.xIonDense[ipHYDROGEN][0],
919 dense.xIonDense[ipHYDROGEN][1],
920 hmi.H2_total,
921 dense.eden,
922 thermal.htot,
923 secondaries.csupra[ipHYDROGEN][0] ,
924 TorF(conv.lgConvIoniz()) ,
925 conv.chConvIoniz());
926 }
927
928 /* this counts number of times we are called by ConvPresTempEdenIoniz,
929 * number of calls in this zone so first call is zero
930 * reset to zero each time ConvPresTempEdenIoniz is called */
931 ++conv.nPres2Ioniz;
932
933 /* this is abort option set with SET PRES IONIZ command,
934 * test on nzone since init can take many iterations
935 * this is seldom used except in special cases */
936 if( conv.nPres2Ioniz > conv.limPres2Ioniz && nzone > 0)
937 {
938 fprintf(ioQQQ," PROBLEM ConvBase sets lgAbort since nPres2Ioniz exceeds limPres2Ioniz. ");
939 fprintf(ioQQQ,"Their values are %li and %li.\n",conv.nPres2Ioniz , conv.limPres2Ioniz);
940 fprintf(ioQQQ," Reset this limit with the SET PRES IONIZ command.\n");
941 lgAbort = true;
942 return 1;
943 }
944
945 /* various checks on the convergence of the current solution */
946 if( eden_sum() )
947 {
948 /* non-zero return indicates abort condition */
949 return 1;
950 }
951
952 /* is electron density converged? */
954 if( fabs(EdenTrue_old - dense.EdenTrue) > fabs(dense.EdenTrue * conv.EdenErrorAllowed/2.) )
955 {
956 conv.setConvIonizFail( "eden chng", EdenTrue_old, dense.EdenTrue);
957 }
958
959 /* check on molecular electron den */
960 if( fabs(EdenFromMolecOld - mole.elec) > fabs(dense.EdenTrue * conv.EdenErrorAllowed/2.) )
961 {
962 conv.setConvIonizFail( "edn chnCO", EdenFromMolecOld, dense.EdenTrue);
963 }
964
965 if( gv.lgGrainElectrons )
966 {
967 /* check on grain electron den */
968 if( fabs(EdenFromGrainsOld - gv.TotalEden) > fabs(dense.EdenTrue * conv.EdenErrorAllowed/2.) )
969 {
970 conv.setConvIonizFail( "edn grn e", EdenFromGrainsOld, gv.TotalEden);
971 }
972
973 /* check on sum of grain and molecular electron den - often two large numbers that cancel */
974 if( fabs( (EdenFromMolecOld-EdenFromGrainsOld) - (mole.elec-gv.TotalEden) ) >
975 fabs(dense.EdenTrue * conv.EdenErrorAllowed/4.) )
976 {
977 conv.setConvIonizFail( "edn mole-grn",
978 (EdenFromMolecOld-EdenFromGrainsOld),
979 (mole.elec-gv.TotalEden));
980 }
981 }
982
983 /* check on heating and cooling if vary temp model
984 * >>chng 08 jul 01, over the code's entire history it had tested whether
985 * this is a constant temperature simulation and did not do this test if
986 * the thermal solution was not done. There are some cases where we do
987 * want to specify the temperature and then find the heating or cooling -
988 * this is done in calculations of cooling curves for instance. With this
989 * change the heating/cooling are converged even in a constant temperature
990 * sim. this does make CT sims run more slowly but with greater accuracy
991 * if heating or cooling is reported */
992 if( !thermal.lgTemperatureConstant )
993 {
994 if( fabs(HeatingOld - thermal.htot)/thermal.htot > conv.HeatCoolRelErrorAllowed/2. )
995 {
996 conv.setConvIonizFail( "heat chg", HeatingOld, thermal.htot);
997 }
998
999 if( fabs(CoolingOld - thermal.ctot)/thermal.ctot > conv.HeatCoolRelErrorAllowed/2. )
1000 {
1001 conv.setConvIonizFail( "cool chg", CoolingOld, thermal.ctot);
1002 }
1003 }
1004
1005 /* check whether molecular abundances are stable */
1006 for( long i=0; i < mole_global.num_calc; ++i )
1007 {
1008 // done relative to total nuclei density so that we can do pure metal plasmas.
1009 if( fabs(mole.species[i].den-mole_save[i])/dense.xNucleiTotal-1. >
1010 conv.EdenErrorAllowed/2.)
1011 {
1012 char chConvIoniz[INPUT_LINE_LENGTH];
1013 sprintf( chConvIoniz, "ch %-4.4s",mole_global.list[i]->label.c_str() );
1014 conv.setConvIonizFail( chConvIoniz,
1015 mole_save[i]/dense.xNucleiTotal,
1016 mole.species[i].den/dense.xNucleiTotal);
1017 }
1018 }
1019
1020 /* >>chng 05 mar 26, add this convergence test - important for molecular or advective
1021 * sims since iso ion solver must sync up with chemistry */
1022 /* remember current ground state population - will check if converged */
1023 for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
1024 {
1025 for( long nelem=ipISO; nelem<LIMELM;++nelem )
1026 {
1027 if( dense.lgElmtOn[nelem] )
1028 {
1029 /* only do this check for "significant" levels of ionization */
1030 /*lint -e644 var possibly not init */
1031 if( dense.xIonDense[nelem][nelem-ipISO]/dense.gas_phase[nelem] > 1e-5 )
1032 {
1033 if( fabs(iso_sp[ipISO][nelem].st[0].Pop()-save_iso_grnd[ipISO][nelem])/SDIV(iso_sp[ipISO][nelem].st[0].Pop())-1. >
1034 conv.EdenErrorAllowed)
1035 {
1036 char chConvIoniz[INPUT_LINE_LENGTH];
1037 sprintf( chConvIoniz,"iso %2li %2li",ipISO, nelem );
1038 conv.setConvIonizFail(chConvIoniz,
1039 save_iso_grnd[ipISO][nelem],
1040 iso_sp[ipISO][nelem].st[0].Pop());
1041 }
1042 }
1043 /*lint +e644 var possibly not init */
1044 }
1045 }
1046 }
1047
1048 /* this counts how many times ConvBase has been called in this iteration,
1049 * located here at bottom of routine so that number is false on first
1050 * call, set to 0 in when iteration starts - used to create itr/zn
1051 * number in printout often used to tell whether this is the very
1052 * first attempt at solution in this iteration */
1053 ++conv.nTotalIoniz;
1054
1055 /* first sweep is over if this is not search phase */
1056 if( !conv.lgSearch )
1057 conv.lgFirstSweepThisZone = false;
1058
1059 /* this was set with STOP NTOTALIONIZ command - only a debugging aid
1060 * by default is zero so false */
1061 if(StopCalc.nTotalIonizStop && conv.nTotalIoniz> StopCalc.nTotalIonizStop )
1062 {
1063 /* non-zero return indicates abort condition */
1064 fprintf(ioQQQ , "ABORT flag set since STOP nTotalIoniz was set and reached.\n");
1065 return 1;
1066 }
1067
1068 if( save.lgTraceConvergeBase )
1069 {
1070 if( iteration > 1 && save.lgTraceConvergeBaseHash )
1071 {
1072 static int iter_punch=-1;
1073 if( iteration !=iter_punch )
1074 fprintf( save.ipTraceConvergeBase, "%s\n",save.chHashString );
1075 iter_punch = iteration;
1076 }
1077
1078 fprintf( save.ipTraceConvergeBase,
1079 "%li\t%.4e\t%.4e\t%.4e\n",
1080 nzone , thermal.htot , thermal.ctot , dense.eden );
1081 }
1082
1083 return 0;
1084}
1085
1086void UpdateUTAs( void )
1087{
1088 DEBUG_ENTRY( "UpdateUTAs()" );
1089
1090 /* only reevaluate this on first pass through on this zone */
1091 if( conv.lgFirstSweepThisZone )
1092 {
1093 /* inner shell ionization */
1094 for( long nelem=0; nelem< LIMELM; ++nelem )
1095 {
1096 for( long ion=0; ion<nelem+1; ++ion )
1097 {
1098 ionbal.UTA_ionize_rate[nelem][ion] = 0.;
1099 ionbal.UTA_heat_rate[nelem][ion] = 0.;
1100 }
1101 }
1102 /* inner shell ionization by UTA lines */
1103 /* this flag turned off with no UTA command */
1104 if( ionbal.lgInnerShellLine_on )
1105 {
1106 for( long i=0; i < nUTA; ++i )
1107 {
1108 /* rateone is inverse lifetime of level against autoionization */
1109 double rateone = UTALines[i].Emis().pump() * UTALines[i].Emis().AutoIonizFrac();
1110 ionbal.UTA_ionize_rate[(*UTALines[i].Hi()).nelem()-1][(*UTALines[i].Hi()).IonStg()-1] += rateone;
1111 /* heating rate in erg atom-1 s-1 */
1112 ionbal.UTA_heat_rate[(*UTALines[i].Hi()).nelem()-1][(*UTALines[i].Hi()).IonStg()-1] += rateone*UTALines[i].Coll().heat();
1113 {
1114 /* DEBUG OTS rates - turn this on to debug line, continuum or both rates */
1115 /*@-redef@*/
1116 enum {DEBUG_LOC=false};
1117 /*@+redef@*/
1118 if( DEBUG_LOC /*&& UTALines[i].nelem==ipIRON+1 && (UTALines[i].IonStg==15||UTALines[i].IonStg==14)*/ )
1119 {
1120 fprintf(ioQQQ,"DEBUG UTA %3i %3i %.3f %.2e %.2e %.2e\n",
1121 (*UTALines[i].Hi()).nelem() , (*UTALines[i].Hi()).IonStg() , UTALines[i].WLAng() ,
1122 rateone, UTALines[i].Coll().heat(),
1123 UTALines[i].Coll().heat()*dense.xIonDense[(*UTALines[i].Hi()).nelem()-1][(*UTALines[i].Hi()).IonStg()-1] );
1124 }
1125 }
1126 }
1127 }
1128 }
1129
1130 return;
1131}
1132
1134{
1135 fixit(); // this routine needs to be enabled
1136 return true;
1137 DEBUG_ENTRY( "lgNetEdenSrcSmall()" );
1138
1139 if( conv.lgSearch )
1140 return true;
1141 fixit(); // grain rates not well tested below
1142 if( gv.lgDustOn() )
1143 return true;
1144
1145 // Check for consistency of explicit source and sink rates for
1146 // electrons with population derived from neutrality
1147 double ionsrc = 0., ionsnk = 0.;
1148 for( long nelem=0; nelem < LIMELM; ++nelem )
1149 {
1150 if( !dense.lgElmtOn[nelem] )
1151 continue;
1152 ionsrc += ionbal.elecsrc[nelem];
1153 ionsnk += ionbal.elecsnk[nelem];
1154 for( long ion_from = 0; ion_from <= nelem + 1; ++ion_from )
1155 {
1156 for( long ion_to = 0; ion_to <= nelem + 1; ++ion_to )
1157 {
1158 if( ion_to-ion_from > 0 )
1159 {
1160 ionsrc += gv.GrainChTrRate[nelem][ion_from][ion_to] *
1161 dense.xIonDense[nelem][ion_from] * (ion_to-ion_from);
1162 }
1163 else if( ion_to-ion_from < 0 )
1164 {
1165 ionsnk += gv.GrainChTrRate[nelem][ion_from][ion_to] *
1166 dense.xIonDense[nelem][ion_from] * (ion_from-ion_to);
1167 }
1168 }
1169 }
1170 }
1171 long ipMElec = findspecies("e-")->index;
1172 const double totsrc = ionsrc + mole.species[ipMElec].src;
1173 const double totsnk = ionsnk + mole.species[ipMElec].snk*dense.EdenTrue;
1174 const double diff = (totsrc - totsnk);
1175 const double ave = ( fabs(totsrc) + fabs(totsnk) )/2.;
1176 fixit(); // Need to tighten up e- population convergence criterion
1177 const double error_allowed = 0.05 * ave; //conv.EdenErrorAllowed * ave;
1178 if( fabs(diff) > error_allowed )
1179 {
1180 enum {DEBUG_LOC=false};
1181 if( DEBUG_LOC )
1182 {
1183 fprintf(ioQQQ,"PROBLEM large NetEdenSrc nzone %li\t%e\t%e\t%e\t%e\n",
1184 nzone,
1185 totsrc/SDIV(totsnk),
1186 dense.EdenTrue,
1187 ionsrc + mole.species[ipMElec].src,
1188 ionsnk + mole.species[ipMElec].snk*dense.EdenTrue);
1189 }
1190 conv.setConvIonizFail( "NetEdenSrc", diff, error_allowed);
1191 return false;
1192 }
1193 else
1194 return true;
1195}
t_atmdat atmdat
Definition atmdat.cpp:6
double ChargTranSumHeat(void)
void ChargTranEval(void)
void atmdat_3body(void)
static double x2[63]
static double x0[83]
static double x1[83]
long int nzone
Definition cddefines.cpp:14
FILE * ioQQQ
Definition cddefines.cpp:7
bool lgAbort
Definition cddefines.cpp:10
long int iteration
Definition cddefines.cpp:16
double fnzone
Definition cddefines.cpp:15
#define ASSERT(exp)
Definition cddefines.h:578
#define STATIC
Definition cddefines.h:97
#define MIN2
Definition cddefines.h:761
const int ipOXYGEN
Definition cddefines.h:312
const int LIMELM
Definition cddefines.h:258
const int INPUT_LINE_LENGTH
Definition cddefines.h:254
char TorF(bool l)
Definition cddefines.h:710
const int NISO
Definition cddefines.h:261
const int ipHELIUM
Definition cddefines.h:306
float realnum
Definition cddefines.h:103
#define MAX2
Definition cddefines.h:782
const int ipNICKEL
Definition cddefines.h:332
bool fp_equal(sys_float x, sys_float y, int n=3)
Definition cddefines.h:812
sys_float SDIV(sys_float x)
Definition cddefines.h:952
const int ipHYDROGEN
Definition cddefines.h:305
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
void fixit(void)
Definition service.cpp:991
int index
Definition mole.h:169
double den
Definition mole.h:358
t_conv conv
Definition conv.cpp:5
@ CONV_BASE_LOOPS
Definition conv.h:76
@ CONV_BASE_ACCELS
Definition conv.h:77
@ CONV_BASE_CALLS
Definition conv.h:75
int eden_sum(void)
Definition eden_sum.cpp:18
int ConvBase(long loopi)
void UpdateUTAs(void)
STATIC bool lgNetEdenSrcSmall(void)
void CoolEvaluate(double *tot)
Definition cool_eval.cpp:45
void HeatZero(void)
Definition heat_sum.cpp:928
void HeatSum(void)
Definition heat_sum.cpp:33
const realnum SMALLFLOAT
Definition cpu.h:191
bool lgElemsConserved(void)
Definition dense.cpp:99
t_dense dense
Definition dense.cpp:24
void lgStatesConserved(long nelem, long ionStage, qList states, long numStates, realnum err_tol, long loop_ion)
Definition dense.cpp:160
t_deuterium deut
Definition deuterium.cpp:8
void SetDeuteriumIonization(const double &xNeutral, const double &xIonized)
Definition deuterium.cpp:26
t_dynamics dynamics
Definition dynamics.cpp:44
void DynaIonize(void)
Definition dynamics.cpp:186
t_elementnames elementnames
void GrainDrive(void)
Definition grains.cpp:1591
GrainVar gv
Definition grainvar.cpp:5
vector< diatomics * > diatoms
Definition h2.cpp:8
vector< diatomics * >::iterator diatom_iter
Definition h2.h:13
void highen(void)
Definition highen.cpp:21
t_hmi hmi
Definition hmi.cpp:5
void ion_recom_calculate(void)
void ion_wrapper(long nelem)
void ion_trim(long int nelem)
Definition ion_trim.cpp:21
t_ionbal ionbal
Definition ionbal.cpp:5
t_iso_sp iso_sp[NISO][LIMELM]
Definition iso.cpp:8
void iso_update_rates(void)
Definition iso_solve.cpp:51
void iso_collapsed_update(void)
Definition iso_solve.cpp:27
void iso_renorm(long nelem, long ipISO, double &renorm)
const int ipH_LIKE
Definition iso.h:62
t_mole_global mole_global
Definition mole.cpp:6
t_mole_local mole
Definition mole.cpp:7
molezone * findspecieslocal(const char buf[])
void mole_drive(void)
void mole_update_sources(void)
molecule * findspecies(const char buf[])
t_opac opac
Definition opacity.cpp:5
void OpacityAddTotal(void)
t_phycon phycon
Definition phycon.cpp:6
void PresTotCurrent(void)
t_rfield rfield
Definition rfield.cpp:8
void RT_line_all(void)
void RT_OTS_PrtRate(double weak, int chFlag)
Definition rt_ots.cpp:712
void RT_OTS(void)
Definition rt_ots.cpp:39
void RT_OTS_Update(double *SumOTS)
Definition rt_ots.cpp:488
t_save save
Definition save.cpp:5
t_secondaries secondaries
t_StopCalc StopCalc
Definition stopcalc.cpp:5
long int nPres2Ioniz
Definition conv.h:152
void setConvIonizFail(const char *reason, double oldval, double newval)
Definition conv.h:107
bool lgElmtOn[LIMELM]
Definition dense.h:146
double xIonDense[LIMELM][LIMELM+1]
Definition dense.h:125
realnum gas_phase[LIMELM]
Definition dense.h:71
char chElementSym[LIMELM][CHARS_ELEMENT_SYM]
TransitionList UTALines("UTALines", &AnonStates)
long int nUTA
Definition taulines.cpp:26
TransitionList TauLines("TauLines", &AnonStates)
t_thermal thermal
Definition thermal.cpp:5
t_trace trace
Definition trace.cpp:5
void DumpLine(const TransitionProxy &t)