cloudy trunk
Loading...
Searching...
No Matches
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-2013 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 convergence failure */
6#include "cddefines.h"
7
8#include "phycon.h"
9#include "rt.h"
10#include "dense.h"
11#include "pressure.h"
12#include "trace.h"
13#include "conv.h"
14#include "pressure_change.h"
15#include "thermal.h"
16#include "geometry.h"
17#include "grainvar.h"
18#include "grains.h"
19
20/*ConvPresTempEdenIoniz solve for current pressure, calls PressureChange, ConvTempEdenIoniz,
21 * called by cloudy
22 * returns 0 if ok, 1 if disaster */
24{
25 DEBUG_ENTRY( "ConvPresTempEdenIoniz()" );
26
27 /* this will count number of times we call ConvBase in this zone,
28 * counter is incremented there
29 * zero indicates first pass through solvers on this zone */
30 conv.nPres2Ioniz = 0;
31 conv.lgFirstSweepThisZone = true;
32 conv.lgLastSweepThisZone = false;
33 /* this will save the history of density - pressure relationship
34 * for the current zone */
35 if( nzone != conv.hist_pres_nzone )
36 {
37 /* first time in this zone - reset history */
38 conv.hist_pres_nzone = nzone;
39 conv.hist_pres_density.clear();
40 conv.hist_pres_current.clear();
41 conv.hist_pres_error.clear();
42 }
43
44 /* should still be true at end */
45 conv.lgConvPops = true;
46
47 /*rel_slope = 0.;*/
48
49 if( trace.nTrConvg>=1 )
50 {
51 fprintf( ioQQQ,
52 " ConvPresTempEdenIoniz1 entered, will call ConvIoniz to initialize\n");
53 }
54
55 /* converge the ionization first, so that we know where we are, and have
56 * a valid foundation to begin the search */
57 /* the true electron density dense.EdenTrue is set in eden_sum called by ConvBase */
58
59 /* this evaluates current pressure, and returns whether or not
60 * it is within tolerance of correct pressure */
61 conv.lgConvPres = false;
62
63 /* convergence trace at this level */
64 if( trace.nTrConvg>=1 )
65 {
66 fprintf( ioQQQ,
67 " ConvPresTempEdenIoniz1 ConvIoniz found following converged: Pres;%c, Temp;%c, Eden;%c, Ion:%c, Pops:%c\n",
68 TorF(conv.lgConvPres),
69 TorF(conv.lgConvTemp),
70 TorF(conv.lgConvEden),
71 TorF(conv.lgConvIoniz()),
72 TorF(conv.lgConvPops));
73 }
74
75 /* trace convergence print at this level */
76 if( trace.nTrConvg>=1 )
77 {
78 fprintf( ioQQQ,
79 "\n ConvPresTempEdenIoniz1 entering main pressure loop.\n");
80 }
81
83
84 if ( strcmp(dense.chDenseLaw,"CPRE") == 0 ||
85 strcmp(dense.chDenseLaw,"DYNA") == 0 )
86 {
87 /* set the initial temperature to the current value, so we will know
88 * if we are trying to jump over a thermal front */
89 double TemperatureInitial = phycon.te;
90
91 /* chng 02 dec 11 rjrw -- ConvIoniz => ConvTempEdenIoniz() here for consistency with inside loop */
92 /* ConvIoniz; */
93 if( ConvTempEdenIoniz() )
94 {
95 return 1;
96 }
97
98 PresMode presmode;
99 presmode.set();
100 solverState st;
101 st.press = pressureZone(presmode);
102
103 /* this will be flag to check for pressure oscillations */
104 bool lgPresOscil = false;
105 /* this is loop where it happened */
106 long nloop_pres_oscil = 0;
107 /* we will use these to check whether hden oscillating - would need to decrease step size */
108 double hden_old = dense.gas_phase[ipHYDROGEN];
109 double hden_chng = 0.;
110 /* this is dP_chng_factor, cut in half when pressure error changes sign */
111 double dP_chng_factor = 1.;
112
113 /* the limit to the number of loops */
114 const int LOOPMAX = 100;
115 /* this will be the limit, which we will increase if no oscillations occur */
116 long LoopMax = LOOPMAX;
117 long loop = 0;
118
119 /* if start of calculation and we are solving for set pressure,
120 * allow a lot more iterations */
121 if( nzone <= 1 && pressure.lgPressureInitialSpecified )
122 LoopMax = 10*LOOPMAX;
123
124 while( (loop < LoopMax) && !conv.lgConvPres && !lgAbort )
125 {
126 /* there can be a pressure or density oscillation early in the search - if not persistent
127 * ok to clear flag */
128 /* >>chng 01 aug 24, if large change in temperature allow lots more loops */
129 if( fabs( TemperatureInitial - phycon.te )/phycon.te > 0.3 )
130 LoopMax = 2*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
145 conv.incrementCounter(PRES_CHANGES);
146 bool lgChange = PressureChange(dP_chng_factor, presmode, st);
147 if( lgChange )
148 {
149 /* heating cooling balance while doing ionization,
150 * this is where the heavy lifting is done, this calls PresTotCurrent,
151 * which sets pressure.PresTotlCurr */
152 int lgTempStatus = ConvTempEdenIoniz();
153 if( lgTempStatus != 0 )
154 {
155 return 1;
156 }
157 }
158
159 /* if product of these two is negative then hden is oscillating */
160 double hden_chng_old = hden_chng;
161 hden_chng = dense.gas_phase[ipHYDROGEN] - hden_old;
162 if( fabs(hden_chng) < SMALLFLOAT )
163 hden_chng = sign( (double)SMALLFLOAT, hden_chng );
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\n",
172 nzone,
173 pressure.PresTotlCurr,
174 pressure.PresTotlError*100.
175 );
176 }
177 }
178
179 /* check whether pressure is oscillating */
180 if( ( ( hden_chng*hden_chng_old < 0. ) ) && loop > 1 )
181 {
182 /*fprintf(ioQQQ,"DEBUG\t%.2f\t%.2e\t%.2e\t%.2e\t%.2e\n",
183 fnzone,pres_chng,pres_chng_old ,hden_chng,hden_chng_old);*/
184 lgPresOscil = true;
185 nloop_pres_oscil = loop;
186 /* >>chng 04 dec 09, go to this factor becoming smaller every time oscillation occurs */
187 dP_chng_factor = MAX2(0.1 , dP_chng_factor * 0.75 );
188 /* dP_chng_factor is how pressure changes with density - pass this to
189 * changing routine if it is stable */
190
191 /*fprintf(ioQQQ,"oscilll %li %.2e %.2e %.2e %.2e dP_chng_factor %.2e\n",
192 loop ,
193 pres_chng,
194 pres_chng_old,
195 hden_chng ,
196 hden_chng_old ,
197 rel_slope);*/
198 }
199
200 /* convergence trace at this level */
201 if( trace.nTrConvg>=1 )
202 {
203 fprintf( ioQQQ,
204 " ConvPresTempEdenIoniz1 %.2f l:%3li nH:%.4e ne:%.4e PCurnt:%.4e err:%6.3f%% dP/dn:%.2e Te:%.4e Osc:%c\n",
205 fnzone,
206 loop,
207 dense.gas_phase[ipHYDROGEN],
208 dense.eden,
209 pressure.PresTotlCurr,
210 /* this is percentage error */
211 100.*pressure.PresTotlError,
212 dP_chng_factor ,
213 phycon.te,
214 TorF(lgPresOscil) );
215 }
216
217 /* increment loop counter */
218 ++loop;
219
220 /* there can be a pressure or density oscillation early in the search - if not persistent
221 * ok to clear flag */
222 if( loop - nloop_pres_oscil > 4 )
223 lgPresOscil = false;
224
225 /* if we hit limit of loop, but no oscillations have happened, then we are
226 * making progress, and can keep going */
227 if( loop == LoopMax && !lgPresOscil )
228 {
229 LoopMax = MIN2( LOOPMAX, LoopMax*2 );
230 }
231 }
232 }
233 else
234 {
235 double targetDensity = zoneDensity();
236 double startingDensity = scalingDensity();
237 double pdelta = conv.MaxFractionalDensityStepPerIteration;
238
239 double logRatio = log(targetDensity/startingDensity);
240 long nstep = (long) ceil(fabs(logRatio)/pdelta);
241 // Ensure at least one update per zone
242 if (nstep == 0)
243 nstep = 1;
244 double density_change_factor = exp(logRatio/nstep);
245
246 for (long i=0; i<nstep; i++)
247 {
248 if (i == nstep - 1)
249 {
250 // Try to ensure exact answer at last iteration
251 density_change_factor = targetDensity/scalingDensity();
252 }
253
254 /* this will evaluate current pressure, update the densities,
255 * determine the wind velocity, and set conv.lgConvPres,
256 * return value is true if density was changed, false if no changes were necessary
257 * if density changes then we must redo the temperature and ionization */
258
259 /* >>chng 04 feb 11, add option to remember current density and pressure */
260 conv.hist_pres_density.push_back( dense.gas_phase[ipHYDROGEN] );
261 conv.hist_pres_current.push_back( pressure.PresTotlCurr );
262 conv.hist_pres_error.push_back( pressure.PresTotlError );
263
264 if( trace.lgTrace )
265 {
266 fprintf( ioQQQ,
267 " DensityUpdate called, changing HDEN from %10.3e to %10.3e Set fill fac to %10.3e\n",
268 dense.gas_phase[ipHYDROGEN], density_change_factor*dense.gas_phase[ipHYDROGEN],
269 geometry.FillFac );
270 }
271
272 ScaleAllDensities((realnum) density_change_factor);
273
274 /* must call TempChange since ionization has changed, there are some
275 * terms that affect collision rates (H0 term in electron collision) */
276 TempChange(phycon.te , false);
277
278 /* heating cooling balance while doing ionization, this is
279 * where the heavy lifting is done, this calls
280 * PresTotCurrent, which sets pressure.PresTotlCurr */
281 if( ConvTempEdenIoniz() )
282 {
283 return 1;
284 }
285
286 /* convergence trace at this level */
287 if( trace.nTrConvg>=1 )
288 {
289 fprintf( ioQQQ,
290 " ConvPresTempEdenIoniz1 %.2f l:%3li nH:%.4e ne:%.4e PCurnt:%.4e err:%6.3f%% Te:%.4e Osc:%c\n",
291 fnzone,
292 i,
293 dense.gas_phase[ipHYDROGEN],
294 dense.eden,
295 pressure.PresTotlCurr,
296 /* this is percentage error */
297 100.*pressure.PresTotlError,
298 phycon.te,
299 TorF(false) );
300 }
301 }
302
304
305 conv.lgConvPres = true;
306 }
307 /* keep track of minimum and maximum temperature */
308 thermal.thist = max((realnum)phycon.te,thermal.thist);
309 thermal.tlowst = min((realnum)phycon.te,thermal.tlowst);
310
311 /* >>chng 04 jan 31, now that all of the physics is converged, determine grain drift velocity */
312 if( gv.lgDustOn() && gv.lgGrainPhysicsOn )
313 GrainDrift();
314
315 /* >>chng 01 mar 14, all ConvFail one more time, no matter how
316 * many failures occurred below. Had been series of if, so multiple
317 * calls per failure possible. */
318 /* >>chng 04 au 07, only announce pres fail here,
319 * we did not converge the pressure */
320 if( !conv.lgConvIoniz() )
321 ConvFail("ioni","");
322 else if( !conv.lgConvEden )
323 ConvFail("eden","");
324 else if( !conv.lgConvTemp )
325 ConvFail("temp","");
326 else if( !conv.lgConvPres )
327 ConvFail("pres","");
328
329 /* this is only a sanity check that the summed continua accurately
330 * reflect all of the individual components. Only include this
331 * when NDEBUG is not set, we are in not debug compile */
332# if !defined(NDEBUG)
333 RT_OTS_ChkSum(0);
334# endif
335
336 return 0;
337}
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 MIN2
Definition cddefines.h:761
T sign(T x, T y)
Definition cddefines.h:800
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
const int ipHYDROGEN
Definition cddefines.h:305
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
long min(int a, long b)
Definition cddefines.h:723
t_conv conv
Definition conv.cpp:5
int ConvTempEdenIoniz(void)
@ PRES_CHANGES
Definition conv.h:82
void ConvFail(const char chMode[], const char chDetail[])
Definition conv_fail.cpp:18
int ConvPresTempEdenIoniz(void)
const realnum SMALLFLOAT
Definition cpu.h:191
void ScaleAllDensities(realnum factor)
Definition dense.cpp:37
t_dense dense
Definition dense.cpp:24
realnum scalingDensity(void)
Definition dense.cpp:378
bool AbundChange()
Definition dense.cpp:269
t_geometry geometry
Definition geometry.cpp:5
void GrainDrift(void)
Definition grains.cpp:4884
GrainVar gv
Definition grainvar.cpp:5
t_phycon phycon
Definition phycon.cpp:6
t_pressure pressure
Definition pressure.cpp:5
void PresTotCurrent(void)
double zoneDensity()
double pressureZone(const PresMode &presmode)
bool PressureChange(double dP_chng_factor, const PresMode &presmode, solverState &st)
void RT_OTS_ChkSum(long int ipPnt)
Definition rt_ots.cpp:631
void TempChange(double TempNew, bool lgForceUpdate)
t_thermal thermal
Definition thermal.cpp:5
t_trace trace
Definition trace.cpp:5