cloudy trunk
Loading...
Searching...
No Matches
conv_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/*ConvEdenIoniz called by ConvTempIonz, calls ConvIoniz solving for eden */
4/*lgConvEden returns true if electron density is converged */
5/*EdenError evaluate ConvIoniz() until ionization has converged and return error on eden */
6#include "cddefines.h"
7#include "dense.h"
8#include "trace.h"
9#include "thermal.h"
10#include "thirdparty.h"
11#include "phycon.h"
12#include "conv.h"
13
14/*lgConvEden returns true if electron density is converged */
15STATIC bool lgConvEden();
16/*EdenError evaluate ConvIoniz() until ionization has converged and return error on eden */
17STATIC double EdenError(double eden);
18
19/*ConvEdenIoniz called by ConvTempEdenIoniz, calls ConvIoniz solving for eden
20 * returns 0 if ok, 1 if abort */
22{
23 DEBUG_ENTRY( "ConvEdenIoniz()" );
24
25 /* this routine is called by ConvTempEdenIoniz, it calls ConvIoniz
26 * and changes the electron density until it converges */
27
28 if( trace.lgTrace )
29 {
30 fprintf( ioQQQ, "\n" );
31 fprintf( ioQQQ, " ConvEdenIoniz entered\n" );
32 }
33 if( trace.nTrConvg>=3 )
34 {
35 fprintf( ioQQQ,
36 " ConvEdenIoniz called, entering eden loop using solver %s.\n",
37 conv.chSolverEden);
38 }
39
40 /* save entry value of eden */
41 double EdenEntry = dense.eden;
42
43 // this branch uses the van Wijngaarden-Dekker-Brent method
44 if( strcmp( conv.chSolverEden , "vWDB" )== 0 )
45 {
46 conv.lgConvEden = false;
47
48 iter_track NeTrack;
49 double n1, error1, n2, error2;
50 // this is the maximum relative step in eden
51 double factor = 0.02;
52 bool lgHysteresis = false;
53
54 for( int n=0; n < 3; ++n )
55 {
56 const int DEF_ITER = 10;
57 // if hysteresis is detected, we should lower the maximum
58 // step to avoid upsetting the lower solvers by large steps
59 if( lgHysteresis )
60 factor /= 5.;
61
62 NeTrack.clear();
63
64 // when dense.EdenTrue becomes negative, error1 > n1 (since n1 > 0.)
65 // a straight copy EdenTrue -> eden would then imply a step > 100% down
66 // all of the code below will cap that to n1*(1.-factor), and eden stays > 0.
67 // this is also asserted in EdenError, the ONLY place where dense.eden is set
68
69 n1 = dense.eden;
70 error1 = EdenError( n1 );
71 NeTrack.add( n1, error1 );
72
73 if( conv.lgSearch && dense.EdenTrue > SMALLFLOAT )
74 n2 = sqrt(dense.eden*dense.EdenTrue);
75 else if( abs(safe_div( error1, n1 )) < factor )
76 n2 = dense.EdenTrue;
77 else
78 n2 = ( error1 > 0. ) ? n1*(1.-factor) : n1*(1.+factor);
79
80 // n1 == n2 will occur if SET EDEN command was given
81 if( !fp_equal( n1, n2 ) )
82 error2 = EdenError( n2 );
83 else
84 error2 = error1;
85 NeTrack.add( n2, error2 );
86
87 int j = 0;
88
89 // now hunt until we have bracketed the solution
90 while( error1*error2 > 0. && j++ < DEF_ITER )
91 {
92 n1 = n2;
93 error1 = error2;
94 double deriv = NeTrack.deriv(5);
95 // the factor 1.2 creates 20% safety margin
96 double step = safe_div( -1.2*error1, deriv, 0. );
97 step = sign( min( abs(step), factor*n1 ), step );
98 n2 = n1 + step;
99 error2 = EdenError( n2 );
100 NeTrack.add( n2, error2 );
101 }
102
103 if( error1*error2 > 0. && trace.nTrConvg >= 3 )
104 {
105 fprintf( ioQQQ, " ConvEdenIoniz: bracket failure 1 n1: %e %e n2: %e %e\n",
106 n1, error1, n2, error2 );
107 NeTrack.print_history();
108 }
109
110 // using the derivative failed, so simply start hunting up or downwards
111 // we may need to take a big step, so max_iter should be big
112 while( error1*error2 > 0. && j++ < 20*DEF_ITER )
113 {
114 n1 = n2;
115 error1 = error2;
116 n2 = ( error1 > 0. ) ? n1*(1.-factor) : n1*(1.+factor);
117 error2 = EdenError( n2 );
118 NeTrack.add( n2, error2 );
119 }
120
121 if( error1*error2 > 0. && trace.nTrConvg >= 3 )
122 {
123 fprintf( ioQQQ, " ConvEdenIoniz: bracket failure 2 n1: %e %e n2: %e %e\n",
124 n1, error1, n2, error2 );
125 NeTrack.print_history();
126 }
127
128 NeTrack.clear();
129
130 // the bracket should have been found, now set up the Brent solver
131 if( NeTrack.init_bracket( n1, error1, n2, error2 ) == 0 )
132 {
133 int nBound = 0;
134
135 // set tolerance to 2 ulp; if bracket gets narrower than 3 ulp we declare
136 // a convergence failure to avoid changes getting lost in machine precision
137 NeTrack.set_tol(2.*DBL_EPSILON*n2);
138
139 double NeNew = 0.5*(n1+n2);
140 for( int i = 0; i < (1<<(n/2))*DEF_ITER; i++ )
141 {
142 // check for convergence, as well as a pathologically narrow bracket
143 if( lgConvEden() || NeTrack.bracket_width() < 3.*DBL_EPSILON*n2 )
144 break;
145
146 NeTrack.add( NeNew, EdenError( NeNew ) );
147 NeNew = NeTrack.next_val(factor);
148
149 // this guards against hysteresis. the symptom of hysteresis is
150 // that EdenTrue ends up being outside the bracket consistently.
151 // if this happens several times in a row, break out of this loop
152 int nVal = NeTrack.in_bounds(dense.EdenTrue);
153 if( nVal == 0 )
154 nBound = 0;
155 else
156 nBound += nVal;
157 if( abs(nBound) >= 3 )
158 {
159 lgHysteresis = true;
160 if( trace.nTrConvg >= 3 )
161 fprintf( ioQQQ, " ConvEdenIoniz: hysteresis detected\n" );
162 break;
163 }
164 }
165 }
166
167 if( conv.lgConvEden )
168 break;
169
170 if( trace.nTrConvg >= 3 )
171 {
172 fprintf( ioQQQ, " ConvEdenIoniz: brent fails\n" );
173 NeTrack.print_history();
174 }
175 }
176
177 if( trace.lgTrace || trace.nTrConvg >= 3 )
178 {
179 fprintf( ioQQQ, " ConvEdenIoniz: entry eden %.4e -> %.4e rel chng %.2f%% accuracy %.2f%%\n",
180 EdenEntry, dense.eden, (safe_div(dense.eden,EdenEntry,1.)-1.)*100.,
181 (safe_div(dense.eden,dense.EdenTrue,1.)-1.)*100. );
182 fprintf( ioQQQ, " ConvEdenIoniz returns converged=%c reason %s\n",
183 TorF(conv.lgConvEden), conv.chConvIoniz() );
184 }
185 }
186 else
187 {
188 fprintf( ioQQQ, "ConvEdenIoniz finds insane solver %s\n", conv.chSolverEden );
189 ShowMe();
190 }
191
192 return 0;
193}
194
195/* returns true if electron density is converged */
197{
198 conv.lgConvEden = ( abs(dense.eden-dense.EdenTrue) < abs(dense.eden)*conv.EdenErrorAllowed );
199 if( !conv.lgConvEden )
200 {
201 conv.setConvIonizFail( "Ne big chg" , dense.EdenTrue, dense.eden);
202 }
203 return conv.lgConvEden;
204}
205
206/* evaluate ConvIoniz() until ionization has converged and return error on eden */
207STATIC double EdenError(double eden)
208{
209 // don't let electron density be zero - logs are taken
210 ASSERT( eden > 0. );
211
212 // this is the only place where the new electron density is set
213 conv.incrementCounter(EDEN_CHANGES);
214 EdenChange( eden );
215
216 for( int i=0; i < 5; ++i )
217 {
218 if( ConvIoniz() )
219 lgAbort = true;
220
221 if( conv.lgConvIoniz() )
222 break;
223 }
224
225 double error = dense.eden - dense.EdenTrue;
226
227 if( trace.nTrConvg >= 3 )
228 fprintf( ioQQQ, " EdenError: eden %.4e EdenTrue %.4e rel. err. %.4e\n",
229 dense.eden, dense.EdenTrue, safe_div(dense.eden,dense.EdenTrue,1.)-1. );
230
231 return error;
232}
233
FILE * ioQQQ
Definition cddefines.cpp:7
bool lgAbort
Definition cddefines.cpp:10
#define ASSERT(exp)
Definition cddefines.h:578
#define STATIC
Definition cddefines.h:97
T sign(T x, T y)
Definition cddefines.h:800
sys_float safe_div(sys_float x, sys_float y, sys_float res_0by0)
Definition cddefines.h:961
char TorF(bool l)
Definition cddefines.h:710
bool fp_equal(sys_float x, sys_float y, int n=3)
Definition cddefines.h:812
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
long min(int a, long b)
Definition cddefines.h:723
void ShowMe(void)
Definition service.cpp:181
int in_bounds(double x) const
Definition iter_track.h:170
void add(double x, double fx)
Definition iter_track.h:120
void print_history() const
Definition iter_track.h:186
double bracket_width() const
Definition iter_track.h:85
void clear()
Definition iter_track.h:76
int init_bracket(double x1, double fx1, double x2, double fx2)
Definition iter_track.h:104
double next_val()
void set_tol(double tol)
Definition iter_track.h:81
double deriv(int n, double &sigma) const
t_conv conv
Definition conv.cpp:5
void EdenChange(double EdenNew)
@ EDEN_CHANGES
Definition conv.h:80
int ConvIoniz(void)
STATIC bool lgConvEden()
STATIC double EdenError(double eden)
int ConvEdenIoniz(void)
const realnum SMALLFLOAT
Definition cpu.h:191
t_dense dense
Definition dense.cpp:24
t_trace trace
Definition trace.cpp:5