cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
conv_pres_temp_eden_ioniz.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and
2  * others. For conditions of distribution and use see copyright notice in license.txt */
3 /*ConvPresTempEdenIoniz solve for current pressure, calls PressureChange, ConvTempEdenIonize,
4  * called by cloudy */
5 /*ConvFail handle conergece failure */
6 #include "cddefines.h"
7 #include "phycon.h"
8 #include "rt.h"
9 #include "dense.h"
10 #include "pressure.h"
11 #include "trace.h"
12 #include "conv.h"
13 #include "grainvar.h"
14 #include "grains.h"
15 
16 /* the limit to the number of loops */
17 /* >>chng 02 jun 13, from 40 to 50 */
18 static const int LOOPMAX = 100;
19 
20 /*ConvPresTempEdenIoniz solve for current pressure, calls PressureChange, ConvTempEdenIoniz,
21  * called by cloudy
22  * returns 0 if ok, 1 if disaster */
24 {
25  long int loop,
26  LoopMax=LOOPMAX;
27  double
28  hden_old ,
29  hden_chng_old ,
30  hden_chng,
31  /*pres_old ,*/
32  /*pres_chng_old ,*/
33  /*pres_chng,*/
34  /*old_slope,*/
35  /*rel_slope,*/
36  dP_chng_factor;
37  bool lgPresOscil;
38  long int nloop_pres_oscil;
39  double TemperatureInitial;
40 
41  DEBUG_ENTRY( "ConvPresTempEdenIoniz()" );
42 
43  /* this will count number of times we call ConvBase in this zone,
44  * counter is incremented there */
45  conv.nPres2Ioniz = 0;
46  loop = 0;
47 
48  /* this will be the limit, which we will increase if no oscillations occur */
49  LoopMax = LOOPMAX;
50  /* set the initial temperature to the current value, so we will know
51  * if we are trying to jump over a thermal front */
52  TemperatureInitial = phycon.te;
53 
54  /* this will be flag to check for pressure oscillations */
55  lgPresOscil = false;
56  /* this is loop where it happened */
57  nloop_pres_oscil = 0;
58  /* should still be true at end */
59  conv.lgConvPops = true;
60 
61  /* we will use these to check whether hden oscillating - would need to decrease step size */
62  hden_old = dense.gas_phase[ipHYDROGEN];
63  hden_chng = 0.;
64  /*pres_old = pressure.PresTotlCurr;*/
65  /*pres_chng = 0.;*/
66 
67  /* this is dP_chng_factor
68  * cut in half when pressure error changes sign */
69  dP_chng_factor = 1.;
70  /*rel_slope = 0.;*/
71 
72  if( trace.nTrConvg>=1 )
73  {
74  fprintf( ioQQQ,
75  " ConvPresTempEdenIoniz1 entered, will call ConvIoniz to initialize\n");
76  }
77 
78  /* converge the ionization first, so that we know where we are, and have
79  * a valid foundation to begin the search */
80  /* the true electron density dense.EdenTrue is set in eden_sum called by ConvBase */
81 
82  /* chng 02 dec 11 rjrw -- ConvIoniz => ConvTempEdenIoniz() here for consistency with inside loop */
83  /* ConvIoniz; */
84  if( ConvTempEdenIoniz() )
85  {
86 
87  return 1;
88  }
89 
90  /* this evaluates current pressure, and returns whether or not
91  * it is within tolerance of correct pressure */
92  conv.lgConvPres = false; /* lgConvPres(); */
93 
94  /* convergence trace at this level */
95  if( trace.nTrConvg>=1 )
96  {
97  fprintf( ioQQQ,
98  " ConvPresTempEdenIoniz1 ConvIoniz found following converged: Pres;%c, Eden;%c, Temp;%c, Ion:%c Pops:%c\n",
99  TorF(conv.lgConvPres) ,
100  TorF(lgConvEden() ),
101  TorF(lgConvTemp() ) ,
103  TorF(conv.lgConvPops));
104  }
105 
106  /* >>chng 01 apr 01, add test for at least 2 loops to get better pressure convergence */
107  /* >>chng 01 oct 31, add test for number of times converged, for constant
108  * pressure will get two valid solutions */
109  /*while( (loop < LoopMax) && !(conv.lgConvPres &&nPresConverged > 1) && !lgAbort )*/
110  /* >>chng 02 dec 12, do not demand two constant pressure being valid - should not be
111  * necessary if first one really is valid, and this caused unneeded second evaluation
112  * in the constant density cases */
113 
114  /* trace convergence print at this level */
115  if( trace.nTrConvg>=1 )
116  {
117  fprintf( ioQQQ,
118  "\n ConvPresTempEdenIoniz1 entering main pressure loop.\n");
119  }
120  while( (loop < LoopMax) && !conv.lgConvPres && !lgAbort )
121  {
122  /* there can be a pressure or density oscillation early in the search - if not persistent
123  * ok to clear flag */
124  /* >>chng 01 aug 24, if large change in temperature allow lots more loops */
125  if( fabs( TemperatureInitial - phycon.te )/phycon.te > 0.3 )
126  LoopMax = 2*LOOPMAX;
127  /* if start of calculation and we are solving for set pressure,
128  * allow a lot more iterations */
130  LoopMax = 10*LOOPMAX;
131 
132  /* change current densities of all constituents if necessary,
133  * PressureChange evaluates lgPresOK, true if pressure is now ok
134  * sets CurrentPressure and CorrectPressure */
135  hden_old = dense.gas_phase[ipHYDROGEN];
136  /*pres_old = pressure.PresTotlCurr;*/
137 
138  /* this will evaluate current pressure, update the densities,
139  * determine the wind velocity, and set conv.lgConvPres,
140  * return value is true if density was changed, false if no changes were necessary
141  * if density changes then we must redo the temperature and ionization
142  * PressureChange contains the logic that determines how to change the density to get
143  * the right pressure */
144  if( PressureChange( dP_chng_factor ) )
145  {
146  /* heating cooling balance while doing ionization,
147  * this is where the heavy lifting is done, this calls PresTotCurrent,
148  * which sets pressure.PresTotlCurr */
149  if( ConvTempEdenIoniz() )
150  {
151  return 1;
152  }
153  }
154 
155  /* if product of these two is negative then hden is oscillating */
156  hden_chng_old = hden_chng;
157  /*pres_chng_old = pres_chng;*/
158  hden_chng = dense.gas_phase[ipHYDROGEN] - hden_old;
159  if( fabs(hden_chng) < SMALLFLOAT )
160  hden_chng = sign( (double)SMALLFLOAT, hden_chng );
161  /*pres_chng = pressure.PresTotlCurr - pres_old;*/
162  /*old_slope = rel_slope;*/
163  /*rel_slope = (pres_chng/pressure.PresTotlCurr) / (hden_chng/dense.gas_phase[ipHYDROGEN]);*/
164 
165  {
166  /*@-redef@*/
167  enum{DEBUG_LOC=false};
168  /*@+redef@*/
169  if( DEBUG_LOC && nzone > 150 && iteration > 1 )
170  {
171  fprintf(ioQQQ,"%li\t%.2e\t%.2e\t%.2e\n",
172  nzone,
176  );
177  }
178  }
179 
180  /* check whether pressure is oscillating */
181  /* >>chng 02 may 31, add check on sign on hden changes */
182  /* >>chng 02 jun 05, add loop > 1 so that don't trigger off old vals that have
183  * not stabilized yet */
184  /* >>chng 04 dec 11, rm test of pressure changing, only check whether hden is changing
185  * very small changes in hden resulted in changes in pressure with some noise due to
186  * finite precision in te solver - that makes very small oscillations that are not
187  * actually oscillations, only noise */
188  if( ( /*( pres_chng*pres_chng_old < 0. )||*/ ( hden_chng*hden_chng_old < 0. ) ) &&
189  loop > 1 )
190  {
191  /*fprintf(ioQQQ,"DEBUG\t%.2f\t%.2e\t%.2e\t%.2e\t%.2e\n",
192  fnzone,pres_chng,pres_chng_old ,hden_chng,hden_chng_old);*/
193  /* the sign of the change in pressure has changed, so things
194  * are oscillating. This would be a problem */
195  lgPresOscil = true;
196  nloop_pres_oscil = loop;
197  /* >>chng 04 dec 09, go to this factor becoming smaller every time oscillation occurs */
198  dP_chng_factor = MAX2(0.1 , dP_chng_factor * 0.75 );
199  /* dP_chng_factor is how pressure changes with density - pass this to
200  * changing routine if it is stable */
201  /* >>chng 04 dec 09, take average of old and new dP_chng_factor to avoid pressure
202  * failures in orion_hii_pdr_pp.in
203  if( loop > 4 && old_slope*rel_slope > 0. )
204  dP_chng_factor = (1.-FRACNEW)*dP_chng_factor + FRACNEW*rel_slope;*/
205 
206  /*fprintf(ioQQQ,"oscilll %li %.2e %.2e %.2e %.2e dP_chng_factor %.2e\n",
207  loop ,
208  pres_chng,
209  pres_chng_old,
210  hden_chng ,
211  hden_chng_old ,
212  rel_slope);*/
213  }
214 
215  /* convergence trace at this level */
216  if( trace.nTrConvg>=1 )
217  {
218  fprintf( ioQQQ,
219  " ConvPresTempEdenIoniz1 %.2f l:%3li nH:%.4e ne:%.4e PCurnt:%.4e PCorct:%.4e err:%6.3f%% dP/dn:%.2e Te:%.4e Osc:%c\n",
220  fnzone,
221  loop,
223  dense.eden,
226  /* this is percentage error */
228  dP_chng_factor ,
229  phycon.te,
230  TorF(lgPresOscil) );
231  }
232 
233  /* increment loop counter */
234  ++loop;
235 
236  /* there can be a pressure or density oscillation early in the search - if not persistent
237  * ok to clear flag
238  * >>chng 04 sep 22, add this logic */
239  if( loop - nloop_pres_oscil > 4 )
240  lgPresOscil = false;
241 
242  /* if we hit limit of loop, but no oscillations have happened, then we are
243  * making progress, and can keep going */
244  if( loop == LoopMax && !lgPresOscil )
245  {
246  LoopMax = MIN2( 100 , LoopMax*2 );
247  }
248  }
249 
250  /* >>chng 04 jan 31, now that all of the physics is converged, determine grain drift velocity */
251  if( gv.lgDustOn && gv.lgGrainPhysicsOn )
252  GrainDrift();
253 
254  /* >>chng 01 mar 14, all ConvFail one more time, no matter how
255  * many failures occurred below. Had been series of if, so multiple
256  * calls per failure possible. */
257  /* >>chng 04 au 07, only announce pres fail here,
258  * we did not converge the pressure */
259  if( !conv.lgConvIoniz )
260  ConvFail("ioni","");
261  else if( !conv.lgConvEden )
262  ConvFail("eden","");
263  else if( !conv.lgConvTemp )
264  ConvFail("temp","");
265  else if( !conv.lgConvPres )
266  ConvFail("pres","");
267 
268  /* this is only a sanity check that the summed continua accurately reflect
269  * all of the individual components. Only include this when NDEBUG is not set,
270  * we are in not debug compile */
271 # if !defined(NDEBUG)
272  RT_OTS_ChkSum(0);
273 # endif
274 
275  return 0;
276 }

Generated for cloudy by doxygen 1.8.4