cloudy trunk
Loading...
Searching...
No Matches
conv_init_solution.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/*ConvInitSolution drive search for initial temperature, for illuminated face */
4/*lgTempTooHigh returns true if temperature is too high */
5/*lgCoolHeatCheckConverge return true if converged, false if change in heating cooling larger than set tolerance */
6#include "cddefines.h"
7#include "physconst.h"
8#include "trace.h"
9#include "struc.h"
10#include "rfield.h"
11#include "mole.h"
12#include "dense.h"
13#include "stopcalc.h"
14#include "heavy.h"
15#include "wind.h"
16#include "geometry.h"
17#include "thermal.h"
18#include "radius.h"
19#include "phycon.h"
20#include "pressure.h"
21#include "conv.h"
22#include "hmi.h"
23#include "dynamics.h"
24
25/* derivative of net cooling wrt temperature to check on sign oscillations */
26static double dCoolNetDTOld = 0;
27
28static double OxyInGrains , FracMoleMax;
29/* determine whether chemistry is important - what is the largest
30 * fraction of an atom in molecules? Also is ice formation
31 * important
32 * sets above two file static variables */
33
34/*lgCoolHeatCheckConverge return true if converged, false if change in heating cooling larger than set tolerance */
36 double *CoolNet,
37 bool lgReset )
38{
39 static double HeatOld=-1. , CoolOld=-1.;
40 DEBUG_ENTRY( "lgCoolHeatCheckConverge()" );
41
42 if( lgReset )
43 {
44 HeatOld = -1;
45 CoolOld = -1;
46 }
47
48 double HeatChange = thermal.htot - HeatOld;
49 double CoolChange = thermal.ctot - CoolOld;
50 /* will use larger of heating or cooling in tests for relative convergence
51 * because in constant temperature models one may have trivially small value */
52 double HeatCoolMax = MAX2( thermal.htot , thermal.ctot );
53
54 /* is the heating / cooling converged?
55 * get max relative change, determines whether converged heating/cooling */
56 double HeatRelChange = fabs(HeatChange)/SDIV( HeatCoolMax );
57 double CoolRelChange = fabs(CoolChange)/SDIV( HeatCoolMax );
58 bool lgConverged = true;
59 if( MAX2( HeatRelChange , CoolRelChange ) > conv.HeatCoolRelErrorAllowed )
60 lgConverged = false;
61
62 *CoolNet = thermal.ctot - thermal.htot;
63
64 HeatOld = thermal.htot;
65 CoolOld = thermal.ctot;
66 return lgConverged;
67}
68
69/* call lgCoolHeatCheckConverge until cooling and heating are converged */
71 double *CoolNet ,
72 double *dCoolNetDT,
73 bool lgReset )
74{
75 const long int LOOP_MAX=20;
76 long int loop = 0;
77 bool lgConvergeCoolHeat = false;
78
79 DEBUG_ENTRY( "lgCoolNetConverge()" );
80
81 while( loop < LOOP_MAX && !lgConvergeCoolHeat && !lgAbort )
82 {
83 if( ConvEdenIoniz() )
84 {
85 /* this is an error return, calculation will immediately stop */
86 lgAbort = true;
87 }
88 lgConvergeCoolHeat = lgCoolHeatCheckConverge( CoolNet, lgReset && loop==0 );
89 if( trace.lgTrace || trace.nTrConvg>=3 )
90 fprintf(ioQQQ," lgCoolNetConverge %li calls to ConvEdenIoniz, converged? %c\n",
91 loop , TorF( lgConvergeCoolHeat ) );
92 ++loop;
93 }
94
95 if( thermal.ConstTemp > 0 )
96 {
97 /* constant temperature model - this is trick so that temperature
98 * updater uses derivative to go to the set constant temperature */
99 *CoolNet = phycon.te - thermal.ConstTemp;
100 *dCoolNetDT = 1.f;
101 }
102 else if( phycon.te < 4000.f )
103 {
104 /* at low temperatures the analytical cooling-heating derivative
105 * is often of no value - the species densities change when the
106 * temperature changes. Use simple dCnet/dT=(C-H)/T - this is
107 * usually too large and results is smaller than optical dT, but
108 * we do eventually converge */
109 *dCoolNetDT = thermal.ctot / phycon.te;
110 }
111 else if( thermal.htot / thermal.ctot > 2.f )
112 {
113 /* we are far from the solution and the temperature is much lower than
114 * equilibrium - analytical derivative from cooling evaluation is probably
115 * wrong since coolants are not correct */
116 *dCoolNetDT = thermal.ctot / phycon.te;
117 }
118 else
119 {
120 *dCoolNetDT = thermal.dCooldT - thermal.dHeatdT;
121 if( dCoolNetDTOld * *dCoolNetDT < 0. )
122 {
123 /* derivative is oscillating */
124 fprintf(ioQQQ,"PROBLEM CoolNet derivative oscillating in initial solution "\
125 "Te=%.3e dCoolNetDT=%.3e CoolNet=%.3e dCooldT=%.3e dHeatdT=%.3e\n",
126 phycon.te , *dCoolNetDT , *CoolNet,
127 thermal.dCooldT , thermal.dHeatdT);
128 *dCoolNetDT = *CoolNet / phycon.te;
129 }
130 }
131 return lgConvergeCoolHeat;
132}
133
134/*ChemImportance find fraction of heavy elements in molecules, O in ices */
136{
137 DEBUG_ENTRY( "ChemImportance()" );
138
139 FracMoleMax = 0.;
140 long int ipMax = -1;
141 for( long i=0; i<mole_global.num_calc; ++i )
142 {
143 if (mole.species[i].location == NULL && !mole_global.list[i]->isMonatomic())
144 {
145 for( molecule::nAtomsMap::iterator atom=mole_global.list[i]->nAtom.begin();
146 atom != mole_global.list[i]->nAtom.end(); ++atom)
147 {
148 long nelem = atom->first->el->Z-1;
149 if( dense.lgElmtOn[nelem] )
150 {
151 realnum dens_elemsp = (realnum) mole.species[i].den*atom->second;
152 if ( FracMoleMax*dense.gas_phase[nelem] < dens_elemsp )
153 {
154 FracMoleMax = dens_elemsp/dense.gas_phase[nelem];
155 ipMax = i;
156 }
157 }
158 }
159 }
160 }
161
162 /* fraction of all oxygen in ices on grains */
163 OxyInGrains = 0;
165 for(long i=0;i<mole_global.num_calc;++i)
166 {
167 if (! mole_global.list[i]->lgGas_Phase && mole_global.list[i]->parentLabel.empty() &&
168 mole_global.list[i]->nAtom.find(elOxygen) != mole_global.list[i]->nAtom.end() )
169 OxyInGrains += mole.species[i].den*mole_global.list[i]->nAtom[elOxygen];
170 }
171 /* this is now fraction of O in ices */
172 OxyInGrains /= SDIV(dense.gas_phase[ipOXYGEN]);
173
174 {
175 /* option to print out entire matrix */
176 enum {DEBUG_LOC=false};
177 if( DEBUG_LOC )
178 {
179 fprintf(ioQQQ,"DEBUG fraction molecular %.2e species %li O ices %.2e\n" ,
180 FracMoleMax , ipMax ,OxyInGrains );
181 }
182 }
183
184 return;
185}
186
188{
189 double TeChangeFactor;
190
191 DEBUG_ENTRY( "FindTempChangeFactor()" );
192
193 /* find fraction of atoms in molecules - need small changes
194 * in temperature if fully molecular since chemistry
195 * network is complex and large changes in solution would
196 * cause linearization to fail */
197 /* this evaluates the global variables FracMoleMax and
198 * OxyInGrains */
200
201 /* Following are series of chemical
202 * tests that determine the factor by which
203 * which can change the temperature. must be VERY small when
204 * ice formation on grains is important due to exponential sublimation
205 * rate. Also small when molecules are important, to keep chemistry
206 * within limits of linearized solver
207 * the complete linearization that is implicit in all these solvers
208 * will not work when large Te jumps occur when molecules/ices
209 * are important - solvers can't return to solution if too far away
210 * keep temp steps small when mole/ice is important */
211 if( OxyInGrains > 0.99 )
212 {
213 TeChangeFactor = 0.999;
214 }
215 else if( OxyInGrains > 0.9 )
216 {
217 TeChangeFactor = 0.99;
218 }
219 /* >>chng 97 mar 3, make TeChangeFactor small when close to lower bound */
220 else if( phycon.te < 5. ||
221 /*>>chng 06 jul 30, or if molecules/ices are important */
222 OxyInGrains > 0.1 )
223 {
224 TeChangeFactor = 0.97;
225 }
226 /*>>chng 07 feb 23, add test of chemistry being very important */
227 else if( phycon.te < 20. || FracMoleMax > 0.1 )
228 {
229 /* >>chng 07 feb 24, changed from 0,9 to 0.95 to converge
230 * pdr_coolbb.in test */
231 TeChangeFactor = 0.95;
232 }
233 else
234 {
235 /* this is the default largest step */
236 TeChangeFactor = 0.8;
237 }
238 return TeChangeFactor;
239}
240
241/*ConvInitSolution drive search for initial temperature, for illuminated face */
243{
244 long int i,
245 ionlup,
246 nelem ,
247 nflux_old,
248 nelem_reflux ,
249 ion_reflux;
250 /* current value of Cooling - Heating */
251 bool lgConvergeCoolHeat;
252
253 double
254 TeChangeFactor,
255 TeBoundLow,
256 TeBoundHigh;
257
258 DEBUG_ENTRY( "ConvInitSolution()" );
259
260 /* set NaN for safety */
261 set_NaN( TeBoundLow );
262 set_NaN( TeBoundHigh );
263
264 /* this counts number of times ConvBase is called by PressureChange, in current zone */
265 conv.nPres2Ioniz = 0;
266 /* this counts how many times ConvBase is called in this iteration
267 * and is flag used by various routines to understand they are
268 * being called the first time*/
269 conv.nTotalIoniz = 0;
270 conv.hist_pres_nzone = -1;
271 conv.hist_temp_nzone = -1;
272 /* ots rates not oscillating */
273 conv.lgOscilOTS = false;
274
275 lgAbort = false;
276 dense.lgEdenBad = false;
277 dense.nzEdenBad = 0;
278 /* these are variables to remember the biggest error in the
279 * electron density, and the zone in which it occurred */
280 conv.BigEdenError = 0.;
281 conv.AverEdenError = 0.;
282 conv.BigHeatCoolError = 0.;
283 conv.AverHeatCoolError = 0.;
284 conv.BigPressError = 0.;
285 conv.AverPressError = 0.;
286
287 /* these are equal if set dr was used, and we know the first dr */
288 if( fp_equal( radius.sdrmin, radius.sdrmax ) )
289 {
290 // lgSdrmaxRel true if sdrmax is relative to current radius, false if limit in cm
291 double rfac = radius.lgSdrmaxRel ? radius.Radius : 1.;
292 radius.drad = MIN2( rfac*radius.sdrmax, radius.drad );
293 radius.drad_mid_zone = radius.drad/2.;
294 radius.drad_x_fillfac = radius.drad * geometry.FillFac;
295 }
296
297 /* the H+ density is zero in sims with no H-ionizing radiation and no cosmic rays,
298 * the code does work in this limit. The H0 density can be zero in limit
299 * of very high ionization where H0 underflows to zero. */
300 ASSERT( dense.xIonDense[ipHYDROGEN][0] >=0 && dense.xIonDense[ipHYDROGEN][1]>= 0.);
301
302 if( trace.nTrConvg )
303 fprintf( ioQQQ, "\nConvInitSolution entered \n" );
304
305 /********************************************************************
306 *
307 * this is second or higher iteration, reestablish original temperature
308 *
309 *********************************************************************/
310 if( iteration != 1 )
311 {
312 /* this is second or higher iteration on multi-iteration model */
313 if( trace.lgTrace || trace.nTrConvg )
314 {
315 fprintf( ioQQQ, " ConvInitSolution called, ITER=%2ld resetting Te to %10.2e\n",
316 iteration, struc.testr[0] );
317 }
318
319 if( trace.lgTrace || trace.nTrConvg )
320 {
321 fprintf( ioQQQ, " search set true\n" );
322 }
323
324 /* search phase must be turned on so that variables such as the ots rates,
325 * secondary ionization, and auger yields, can converge more quickly to
326 * proper values */
327 conv.lgSearch = true;
328 conv.lgFirstSweepThisZone = true;
329 conv.lgLastSweepThisZone = false;
330
331 /* reset to the zone one temperature from the previous iteration */
332 TempChange( struc.testr[0] , false);
333
334 /* find current pressure - sets pressure.PresTotlCurr */
336
337 /* get new initial temperature and pressure */
339 {
340 /* this is an error return, calculation will immediately stop */
341 lgAbort = true;
342 return lgAbort;
343 }
344
345 if( trace.lgTrace || trace.nTrConvg )
346 {
347 fprintf( ioQQQ, " ========================================================================================================================\n");
348 fprintf( ioQQQ, " ConvInitSolution: search set false 2 Te=%.3e========================================================================================\n" , phycon.te);
349 fprintf( ioQQQ, " ========================================================================================================================\n");
350 }
351 conv.lgSearch = false;
352
353 /* get temperature a second time */
354 if( ConvTempEdenIoniz() )
355 {
356 /* this is an error return, calculation will immediately stop */
357 lgAbort = true;
358 return lgAbort;
359 }
360
361 /* the initial pressure should now be valid
362 * this sets values of pressure.PresTotlCurr */
364 }
365
366 else
367 {
368 /********************************************************************
369 *
370 * do first te from scratch
371 *
372 *********************************************************************/
373
374 // Set to false to switch to using only standard temperature convergence method
375 const bool lgDoInitConv = true;
376 /* say that we are in search phase */
377 conv.lgSearch = true;
378 conv.lgFirstSweepThisZone = true;
379 conv.lgLastSweepThisZone = false;
380
381 if( trace.lgTrace )
382 {
383 fprintf( ioQQQ, " ConvInitSolution called, new temperature.\n" );
384 }
385
386 /* coming up to here Te is either 4,000 K (usually) or 10^6 */
387 TeBoundLow = phycon.TEMP_LIMIT_LOW/1.001;
388
389 /* set temperature floor option - StopCalc.TeFloor is usually
390 * zero but reset this this command - will go over to constant
391 * temperature calculation if temperature falls below floor */
392 TeBoundLow = MAX2( TeBoundLow , StopCalc.TeFloor );
393
394 /* highest possible temperature - must not get this high since
395 * parts of code will abort if value is reached.
396 * divide by 1.2 to prevent checking on temp > 1e10 */
397 TeBoundHigh = phycon.TEMP_LIMIT_HIGH/1.2;
398
399 /* set initial temperature, options are constant temperature,
400 * or approach equilibrium from either high or low TE */
401 double TeFirst;
402 if( thermal.ConstTemp > 0 )
403 {
404 /* constant temperature , set 4000 K then relax to desired value
405 * for cold temperatures, to slowly approach fully molecular solution.
406 * if constant temperature is higher than 4000., go right to it */
407 /* true allow ionization stage trimming, false blocks it */
408 TeFirst = thermal.ConstTemp ;
409 if (lgDoInitConv)
410 TeFirst = MAX2( 4000. , TeFirst );
411
412 /* override this if molecule deliberately disabled */
413 if( mole_global.lgNoMole )
414 TeFirst = thermal.ConstTemp;
415 }
416 else if( thermal.lgTeHigh )
417 {
418 /* approach from high TE */
419 TeFirst = MIN2( 1e6 , TeBoundHigh );
420 }
421 else
422 {
423 /* approach from low TE */
424 TeFirst = MAX2( 4000., TeBoundLow );
425 }
426
427 // initial kinetic temperature should be at or above the local
428 // energy density temperature if the radiation field is well
429 // coupled to the matter, but never let this overrule
430 // constant temperature command
431 if( !thermal.lgTemperatureConstant )
432 TeFirst = max( TeFirst , phycon.TEnerDen );
433
434 TempChange(TeFirst , false);
435 if( trace.lgTrace || trace.nTrConvg>=2 )
436 fprintf(ioQQQ,"ConvInitSolution doing initial solution with temp=%.2e\n",
437 phycon.te);
438
439 if (lgDoInitConv)
440 {
441 /* this sets values of pressure.PresTotlCurr */
443
444 thermal.htot = 1.;
445 thermal.ctot = 1.;
446
447 /* call cooling, heating, opacity, loop to convergence
448 * this is very first call to it, by default is at 4000K */
449
450 double CoolNet=0, dCoolNetDT=0;
451 /* do ionization twice to get stable solution, evaluating initial dR each time */
452 lgConvergeCoolHeat = false;
453 for( ionlup=0; ionlup<2; ++ionlup )
454 {
455 if( trace.lgTrace || trace.nTrConvg>=2 )
456 fprintf( ioQQQ, "ConvInitSolution calling lgCoolNetConverge "
457 "initial setup %li with Te=%.3e C=%.2e H=%.2e d(C-H)/dT %.3e\n",
458 ionlup,phycon.te , thermal.ctot , thermal.htot,dCoolNetDT );
459 /* do not flag oscillating d(C-H)/dT until stable */
460 dCoolNetDTOld = 0.;
461 lgConvergeCoolHeat = lgCoolNetConverge( &CoolNet , &dCoolNetDT, true );
462 if( lgAbort )
463 break;
464 /* set thickness of first zone, this affects the pumping rates due
465 * to correction for attenuation across zone, so must be known
466 * for ConvEdenIoniz to get valid solution - this is why it
467 * is in this loop */
468 radius_first();
469 }
470 if( !lgConvergeCoolHeat )
471 fprintf(ioQQQ," PROBLEM ConvInitSolution did not converge the heating-cooling during initial search.\n");
472
473 if( lgAbort )
474 {
475 /* we hit an abort */
476 fprintf( ioQQQ, " DISASTER PROBLEM ConvInitSolution: Search for "
477 "initial conditions aborted - lgAbort set true.\n" );
478 ShowMe();
479 /* this is an error return, calculation will immediately stop */
480 return lgAbort;
481 }
482
483 /* we now have error in C-H and its derivative - following is better
484 * derivative for global case where we may be quite far from C==H */
485 if( thermal.ConstTemp<=0 )
486 dCoolNetDT = thermal.ctot / phycon.te;
487 bool lgConvergedLoop = false;
488 long int LoopThermal = 0;
489 /* initial temperature is usually 4000K, if the dTe scale factor is 0.95
490 * it will take 140 integrations to lower temperature to 3K,
491 * and many more if ices are important. */
492 const long int LIMIT_THERMAL_LOOP=300;
493 double CoolMHeatSave=-1. , TempSave=-1. , TeNew=-1.,CoolSave=-1;
494 while( !lgConvergedLoop && LoopThermal < LIMIT_THERMAL_LOOP )
495 {
496 /* change temperature until sign of C-H changes */
497 CoolMHeatSave = CoolNet;
498 dCoolNetDTOld = dCoolNetDT;
499 CoolSave = thermal.ctot;
500 TempSave = phycon.te;
501
502 /* find temperature change factor, a number less than 1*/
503 TeChangeFactor = FindTempChangeFactor();
504 ASSERT( TeChangeFactor>0. && TeChangeFactor< 1. );
505
506 TeNew = phycon.te - CoolNet / dCoolNetDT;
507
508 TeNew = MAX2( phycon.te*TeChangeFactor , TeNew );
509 TeNew = MIN2( phycon.te/TeChangeFactor , TeNew );
510 /* change temperature */
511 TempChange(TeNew , true);
512 lgConvergeCoolHeat = lgCoolNetConverge( &CoolNet , & dCoolNetDT, false );
513
514 if( trace.lgTrace || trace.nTrConvg>=2 )
515 fprintf(ioQQQ,"ConvInitSolution %4li TeNnew=%.3e "
516 "TeChangeFactor=%.3e Cnet/H=%.2e dCnet/dT=%10.2e Ne=%.3e"
517 " Ograins %.2e FracMoleMax %.2e\n",
518 LoopThermal , TeNew , TeChangeFactor ,
519 CoolNet/SDIV(thermal.htot) , dCoolNetDT,
520 dense.eden , OxyInGrains , FracMoleMax );
521
522 if( lgAbort )
523 return lgAbort;
524
525 /* keep changing temperature until sign of CoolNet changes
526 * in constant temperature case CoolNet is delta Temperature
527 * so is zero for last pass in this loop
528 * this is for constant temperature case */
529 if( fabs(CoolNet)< SMALLDOUBLE )
530 /* CoolNet is zero */
531 lgConvergedLoop = true;
532 else if( (CoolNet/fabs(CoolNet))*(CoolMHeatSave/fabs(CoolMHeatSave)) <= 0)
533 /* change in sign of net cooling */
534 lgConvergedLoop = true;
535 else if( phycon.te <= TeBoundLow || phycon.te>=TeBoundHigh)
536 lgConvergedLoop = true;
537 ++LoopThermal;
538 }
539
540 if( !lgConvergedLoop )
541 fprintf(ioQQQ,"PROBLEM ConvInitSolution: temperature solution not "
542 "found in initial search, final Te=%.2e\n",
543 phycon.te );
544
545 /* interpolate temperature where C=H if not constant temperature sim
546 * will have set constant temperature mode above if we hit temperature
547 * floor */
548 if( thermal.ConstTemp<=0 &&
549 ! fp_equal( TempSave , phycon.te ) )
550 {
551 if( trace.lgTrace || trace.nTrConvg>=2 )
552 fprintf(ioQQQ," Change in sign of C-H found, Te brackets %.3e "
553 "%.3e Cool brackets %.3e %.3e ",
554 TempSave , phycon.te , CoolMHeatSave , CoolNet);
555 /* interpolate new temperature assuming Cool = a T^alpha,
556 * first find alpha */
557 double alpha = log(thermal.ctot/CoolSave) / log( phycon.te/TempSave);
558 if( fabs(alpha) < SMALLFLOAT )
559 /* alpha close to 0 means constant temperature */
560 TeNew = phycon.te;
561 else
562 {
563 /* next find log a - work in logs to avoid infinities */
564 if( thermal.ctot<0. || thermal.htot<0. )
566 double Alog = log( thermal.ctot ) - alpha * log( phycon.te );
567 /* the interpolated temperature where heating equals cooling */
568 double TeLn = (log( thermal.htot ) - Alog) / alpha ;
569
570 /* a sanity check */
571 if( TeLn < log( MIN2(phycon.te , TempSave )) )
572 TeNew = MIN2(phycon.te , TempSave );
573 else if( TeLn > log( MAX2(phycon.te , TempSave )) )
574 TeNew = MAX2(phycon.te , TempSave );
575 else
576 TeNew = pow( EE , TeLn );
577 }
578
579 ASSERT( TeNew>=MIN2 ( TempSave , phycon.te ) ||
580 TeNew<=MAX2 ( TempSave , phycon.te ));
581
582 if( trace.lgTrace || trace.nTrConvg>=2 )
583 fprintf(ioQQQ," interpolating to Te %.3e \n\n",
584 TeNew);
585
586 /* change temperature */
587 TempChange(TeNew , true);
588 }
589 }
590
591 if( ConvTempEdenIoniz() )
592 {
593 /* this is an error return, calculation will immediately stop */
594 lgAbort = true;
595 return lgAbort;
596 }
597
598 conv.lgSearch = false;
599
600 // Solve again without search phase lo-fi physics before starting on
601 // anything which might have a long-term impact
602 if( ConvTempEdenIoniz() )
603 {
604 /* this is an error return, calculation will immediately stop */
605 lgAbort = true;
606 return lgAbort;
607 }
608
609 /* this sets values of pressure.PresTotlCurr */
611
612 if( trace.lgTrace || trace.nTrConvg )
613 {
614 fprintf( ioQQQ, "\nConvInitSolution: end 1st iteration search phase "
615 "finding Te=%.3e\n\n" , phycon.te);
616 }
617
618 if( trace.lgTrace )
619 {
620 fprintf( ioQQQ, " ConvInitSolution return, TE:%10.2e==================\n",
621 phycon.te );
622 }
623 }
624
625 /* we now have a fairly good temperature and ionization
626 * iteration is 1 on first iteration - remember current pressure
627 * on first iteration so we can hold this constant in constant
628 * pressure simulation.
629 *
630 * flag dense.lgDenseInitConstant false if
631 * *constant pressure reset* is used -
632 * default is true, after first iteration initial density is used for
633 * first zone no matter what pressure results. Small changes occur due
634 * to radiative transfer convergence,
635 * when set false with *reset* option the density on later iterations
636 * can change to keep the pressure constant */
637 if( iteration==1 || dense.lgDenseInitConstant )
638 {
639 double PresNew = pressure.PresTotlCurr;
640
641 /* option to specify the pressure as option on constant pressure command */
642 if( pressure.lgPressureInitialSpecified )
643 /* this is log of nT product - if not present then set zero */
644 PresNew = pressure.PressureInitialSpecified;
645 if( trace.lgTrace )
646 {
647 fprintf( ioQQQ,
648 " PresTotCurrent 1st zn old values of PresTotlInit:%.3e"
649 " to %.3e \n",
650 pressure.PresTotlInit,
651 PresNew);
652 }
653
654 pressure.PresTotlInit = PresNew;
655 }
656
657 if( dynamics.lgTimeDependentStatic )
658 {
659 // some sort of time dependent sim
660 static double PresTotalInitTime;
661 if( iteration <= dynamics.n_initial_relax )
662 {
663 PresTotalInitTime = pressure.PresTotlInit;
664 }
665 else if (dense.lgPressureVaryTime)
666 {
667 pressure.PresTotlInit = PresTotalInitTime *
668 pow( 1.+(dynamics.time_elapsed/dense.PressureVaryTimeTimescale) ,
669 dense.PressureVaryTimeIndex);
670 }
671// fprintf(ioQQQ,"DEBUG conv_init_solution time dependent time=%.2e sets "
672// "pressure to %.2e\n", dynamics.time_elapsed ,
673// pressure.PresTotlInit);
674 }
675
676 /* now bring current pressure to the correct pressure - must do this
677 * before initial values are saved in iter start/end */
679
680 /* this counts number of times ConvBase is called by PressureChange, in
681 * current zone these are reset here, so that we count from first zone
682 * not search */
683 conv.nPres2Ioniz = 0;
684
685 dense.lgEdenBad = false;
686 dense.nzEdenBad = 0;
687
688 /* save counter upon exit so we can see how many iterations were
689 * needed to do true zones */
690 conv.nTotalIoniz_start = conv.nTotalIoniz;
691
692 /* these are variables to remember the biggest error in the
693 * electron density, and the zone in which it occurred */
694 conv.BigEdenError = 0.;
695 conv.AverEdenError = 0.;
696 conv.BigHeatCoolError = 0.;
697 conv.AverHeatCoolError = 0.;
698 conv.BigPressError = 0.;
699 conv.AverPressError = 0.;
700
701 /*>>chng 06 jun 09, only reset on first try - otherwise have stable solution? */
702 if( iteration == 1 )
703 {
704 /* now remember some things we may need even in first zone, these
705 * are normally set towards end of zone calculation in RT_tau_inc */
706 struc.testr[0] = (realnum)phycon.te;
707 /* >>chng 02 May 2001 rjrw: add hden for dilution */
708 struc.hden[0] = dense.gas_phase[ipHYDROGEN];
709 /* pden is the total number of particles per unit vol */
710 struc.DenParticles[0] = dense.pden;
711 struc.heatstr[0] = thermal.htot;
712 struc.coolstr[0] = thermal.ctot;
713 struc.volstr[0] = (realnum)radius.dVeffAper;
714 struc.drad_x_fillfac[0] = (realnum)radius.drad_x_fillfac;
715 struc.histr[0] = dense.xIonDense[ipHYDROGEN][0];
716 struc.hiistr[0] = dense.xIonDense[ipHYDROGEN][1];
717 struc.ednstr[0] = (realnum)dense.eden;
718 struc.o3str[0] = dense.xIonDense[ipOXYGEN][2];
719 struc.DenMass[0] = dense.xMassDensity;
720 struc.drad[0] = (realnum)radius.drad;
721 }
722
723 /* check that nflux extends above IP of highest ionization species present.
724 * for collisional case it is possible for species to exist that are higher
725 * IP than the limit to the continuum. Need continuum to encompass all
726 * possible emission - to account for diffuse emission
727 * NB
728 * on second iteration of multi-iteration model this may result in rfield.nflux increasing
729 * which can change the solution */
730 nflux_old = rfield.nflux;
731 nelem_reflux = -1;
732 ion_reflux = -1;
733 for( nelem=2; nelem < LIMELM; nelem++ )
734 {
735 /* do not include hydrogenic species in following */
736 for( i=0; i < nelem; i++ )
737 {
738 if( dense.xIonDense[nelem][i+1] > 0. )
739 {
740 if( Heavy.ipHeavy[nelem][i] > rfield.nflux )
741 {
742 rfield.nflux = Heavy.ipHeavy[nelem][i];
743 nelem_reflux = nelem;
744 ion_reflux = i;
745 }
746 }
747 }
748 }
749
750 /* was the upper limit to the continuum updated? if so need to define
751 * continuum variables to this new range */
752 if( nflux_old != rfield.nflux )
753 {
754 /* zero out parts of rfield arrays that were previously undefined */
755 rfield_opac_zero( nflux_old-1 , rfield.nflux );
756
757 /* if continuum reset up, we need to define gaunt factors through high end */
758 /*tfidle(false);*/
759 /* this calls tfidle, among other things */
760 /* this sets values of pressure.PresTotlCurr */
762
763 /* redo ionization and update opacities */
764 if( ConvBase(1) )
765 {
766 /* this is catastrophic failure */
767 lgAbort = true;
768 return lgAbort;
769 }
770
771 /* need to update continuum opacities */
772 if( trace.lgTrace )
773 {
774 fprintf(ioQQQ," nflux updated from %li to %li, anu from %g to %g \n",
775 nflux_old , rfield.nflux ,
776 rfield.anu[nflux_old] , rfield.anu[rfield.nflux] );
777 fprintf(ioQQQ," caused by element %li ion %li \n",
778 nelem_reflux ,ion_reflux );
779 }
780 }
781
782 /* dynamics solution - change density slightly
783 * and call pressure solver to see if it returns to original density */
784 if( strcmp(dense.chDenseLaw,"DYNA") == 0 )
785 {
786 /* >>chng 04 apr 23, change pressure and let solver come back to correct
787 * pressure. This trick makes sure we are correctly either super or sub sonic
788 * since solver will jump from one branch to the other if necessary */
789 static const double PCHNG = 0.98;
790 /* this ignores return condition, assume must be ok if got this far */
792 {
793 /* this is an error return, calculation will immediately stop */
794 lgAbort = true;
795 return lgAbort;
796 }
797
798 pressure.PresTotlInit *= PCHNG;
800 {
801 /* this is an error return, calculation will immediately stop */
802 lgAbort = true;
803 return lgAbort;
804 }
805
806 pressure.PresTotlInit /= PCHNG;
808 {
809 /* this is an error return, calculation will immediately stop */
810 lgAbort = true;
811 return lgAbort;
812 }
813
814# undef PCHNG
815 }
816 /* this is the only valid return and lgAbort should be false*/
817 return lgAbort;
818}
FILE * ioQQQ
Definition cddefines.cpp:7
bool lgAbort
Definition cddefines.cpp:10
long int iteration
Definition cddefines.cpp:16
#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
char TorF(bool l)
Definition cddefines.h:710
float realnum
Definition cddefines.h:103
#define MAX2
Definition cddefines.h:782
long max(int a, long b)
Definition cddefines.h:775
bool fp_equal(sys_float x, sys_float y, int n=3)
Definition cddefines.h:812
NORETURN void TotalInsanity(void)
Definition service.cpp:886
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 ShowMe(void)
Definition service.cpp:181
void rfield_opac_zero(long lo, long ihi)
t_conv conv
Definition conv.cpp:5
int ConvBase(long loopi)
int ConvTempEdenIoniz(void)
int ConvPresTempEdenIoniz(void)
int ConvEdenIoniz(void)
static double OxyInGrains
STATIC bool lgCoolHeatCheckConverge(double *CoolNet, bool lgReset)
STATIC bool lgCoolNetConverge(double *CoolNet, double *dCoolNetDT, bool lgReset)
static double dCoolNetDTOld
bool ConvInitSolution()
STATIC void ChemImportance(void)
static double FracMoleMax
double FindTempChangeFactor(void)
void CoolSave(FILE *io, char chJob[])
Definition cool_save.cpp:20
void set_NaN(sys_float &x)
Definition cpu.cpp:682
const double SMALLDOUBLE
Definition cpu.h:195
const realnum SMALLFLOAT
Definition cpu.h:191
t_dense dense
Definition dense.cpp:24
t_dynamics dynamics
Definition dynamics.cpp:44
t_geometry geometry
Definition geometry.cpp:5
static const long LOOP_MAX
t_Heavy Heavy
Definition heavy.cpp:5
t_mole_global mole_global
Definition mole.cpp:6
t_mole_local mole
Definition mole.cpp:7
ChemAtomList unresolved_atom_list
t_phycon phycon
Definition phycon.cpp:6
UNUSED const double EE
Definition physconst.h:23
t_pressure pressure
Definition pressure.cpp:5
void PresTotCurrent(void)
t_radius radius
Definition radius.cpp:5
void radius_first(void)
t_rfield rfield
Definition rfield.cpp:8
t_StopCalc StopCalc
Definition stopcalc.cpp:5
t_struc struc
Definition struc.cpp:6
void TempChange(double TempNew, bool lgForceUpdate)
t_thermal thermal
Definition thermal.cpp:5
t_trace trace
Definition trace.cpp:5