cloudy trunk
Loading...
Searching...
No Matches
newton_step.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/*mole_newton_step line-search improvement in chemical network */
4/* >> chng 02 nov 7 rjrw, Mole Moreliano:
5 * changes to linearized iterative form */
6/*lint -e778 const express evaluates to 0 */
7/*lint -e725 expect positive indentation */
8#include "cddefines.h"
9#include "newton_step.h"
10#include "conv.h"
11#include "thirdparty.h"
12#include "mole.h"
13#include "mole_priv.h"
14
15#define ABSLIM 1e-12
16#define ERRLIM 1e-12
17# ifdef MAT
18# undef MAT
19# endif
20# define MAT(a,I_,J_) ((a)[(I_)*(n)+(J_)])
21
22STATIC void mole_system_error(long n, long merror,
23 const valarray<double> &a,
24 const valarray<double> &b);
25
26enum {PRINTSOL = false};
27
28/* mole_newton_step -- improve balance in chemical network along
29 * descent direction, step limited to ensure improvement */
30bool newton_step(GroupMap &MoleMap, const valarray<double> &b0vec, valarray<double> &b2vec,
31 realnum *eqerror, realnum *error, const long n, double *rlimit, double *rmax,
32 valarray<double> &escale,
33 void (*jacobn)(GroupMap &MoleMap,
34 const valarray<double> &b2vec,
35 double * const ervals, double * const amat,
36 const bool lgJac, bool *lgConserve))
37{
38 bool lgOK=true;
39 long int i,
40 loop;
41
42 double
43 f1,
44 f2,
45 error2,
46 error1,
47 error0,
48 erroreq,
49 grad,
50 pred;
51
52 valarray<double>
53 amat(n*n),
54 b1vec(n),
55 ervals(n),
56 ervals1(n),
57 ervals0(n);
58 int32 merror;
59
60 DEBUG_ENTRY( "newton_step()" );
61
62 for( i=0; i < n; i++ )
63 {
64 b2vec[i] = b0vec[i];
65 ervals0[i] = 0.0;
66 ervals1[i] = 0.0;
67 }
68
69 pred = error0 = error2 = erroreq = grad = f1 = f2 = 0.;
70
71 conv.incrementCounter(NEWTON);
72 const int LOOPMAX = 40;
73 for (loop=0;loop<LOOPMAX;loop++)
74 {
75 bool lgConserve;
76 jacobn(MoleMap, b2vec,get_ptr(ervals),get_ptr(amat),loop==0,&lgConserve);
77
78 if (loop == 0)
79 {
80 for( i=0; i < n; i++ )
81 {
82 escale[i] = MAT(amat,i,i);
83 }
84 // First time through this case, set rlimit by guesswork based on
85 // maintaining matrix condition
86 *rmax = 0.0;
87 for( i=0; i<n; ++i )
88 {
89 if (-escale[i]>*rmax)
90 {
91 *rmax = -escale[i];
92 }
93 }
94 if (*rlimit < 0.0)
95 {
96 // large generally more stable, small more quickly convergent
97 *rlimit = 1e-19 * (*rmax);
98 }
99 else if (*rlimit > *rmax)
100 {
101 // large generally more stable, small more quickly convergent
102 *rlimit = *rmax;
103 }
104 for( i=0; i < n; i++ )
105 {
106 if (! lgConserve || ( (! groupspecies[i]->isMonatomic()) || groupspecies[i]->charge < 0 ) )
107 {
108 // Only apply *rlimit to rows which haven't been
109 // overwritten by conservation constraints
110 MAT(amat,i,i) -= *rlimit;
111 }
112 }
113 }
114
115 // External error limit based on difference to state relaxed to
116 // equilibrium
117 erroreq = 0.;
118 for( i=0; i < n; i++ )
119 {
120 double etmp = ervals[i]/(SMALLABUND*(*rmax)+fabs(b0vec[i]*escale[i]));
121 erroreq += etmp*etmp;
122 }
123
124 // Internal (step) limit based on finite-(pseudo-)timestep solution
125 for (i=0; i < n; i++)
126 {
127 ervals1[i] = ervals[i];
128 if (! lgConserve || ( (! groupspecies[i]->isMonatomic()) || groupspecies[i]->charge < 0 ) )
129 ervals1[i] -= (*rlimit)*(b2vec[i]-b0vec[i]);
130 }
131
132 error1 = 0.;
133 long maxi = -1;
134 double emax0 = 0.0, emax1 = 0.0;
135 for( i=0; i < n; i++ )
136 {
137 if (! lgConserve || ( (! groupspecies[i]->isMonatomic()) || groupspecies[i]->charge < 0 ) )
138 {
139 /* Scale the errors weighting to ensure trace species are accurate */
140 double etmp = ervals1[i]/(SMALLABUND*(*rmax)+fabs(b0vec[i]*escale[i]));
141 double etmp0 = ervals0[i]/(SMALLABUND*(*rmax)+fabs(b0vec[i]*escale[i]));
142 etmp *= etmp;
143 etmp0 *= etmp0;
144 if (fabs(etmp-etmp0) > fabs(emax1-emax0) )
145 {
146 maxi = i;
147 emax1 = etmp;
148 emax0 = etmp0;
149 }
150 error1 += etmp;
151 }
152 }
153
154 if (loop == 0)
155 {
156 for( i=0; i < n; i++ )
157 {
158 ervals0[i] = ervals1[i];
159 b1vec[i] = ervals1[i];
160 }
161 merror = solve_system(amat,b1vec,n,mole_system_error);
162
163 if (merror != 0) {
164 *eqerror = *error = 1e30f;
165 lgOK = false;
166 return lgOK;
167 }
168 error0 = error1;
169 grad = -2*error0;
170 f1 = 1.0;
171
172 } else {
173 //fprintf(ioQQQ,"Newt1 %ld stp %11.4g err %11.4g erreq %11.4g rlimit %11.4g grd %11.4g fdg %11.4g e0 %11.4g de %11.4g\n",
174 // loop,f1,error1,erroreq,*rlimit,grad,(error1-error0)/f1,error0,error1-error0);
175 if (error1 < (1-2e-4*f1)*error0 || error1 < 1e-20)
176 {
177 break;
178 }
179 // Backtrack using quadratic or cubic fit, ref Dennis & Schnabel 1996
180 if (loop == 1)
181 {
182 f2 = f1;
183 f1 *= -0.5*f1*grad/(error1-error0-f1*grad);
184 pred = error0+0.5*grad*f1;
185 }
186 else
187 {
188 if (0)
189 {
190 fprintf(ioQQQ,"Newt %ld stp %11.4g err %11.4g erreq %11.4g rlimit %11.4g grd %11.4g fdg %11.4g e0 %11.4g de %11.4g pred %11.4g\n",
191 loop,f1,error1,erroreq,*rlimit,grad,(error1-error0)/f1,error0,error1-error0,pred
192 );
193 if (maxi != -1)
194 fprintf(ioQQQ,"Maxi %ld %s emax from %11.4g of %11.4g -> %11.4g of %11.4g, diff %11.4g\n",
195 maxi,groupspecies[maxi]->label.c_str(),emax0,error0, emax1,error1,emax1-emax0);
196 }
197 // NB we must be careful how a is calculated because terms can be nearly equal and/or
198 // vastly different from other terms.
199 double a = (error1-error0)/(f1*f1) - (error2-error0)/(f2*f2);
200 a += grad * (1./f2 - 1./f1);
201 // similarly calculate b
202 double b = -f2*(error1-error0)/(f1*f1) + f1*(error2-error0)/(f2*f2);
203 b += grad * (f2/f1 - f1/f2);
204 double ft = f1;
205 //fprintf(ioQQQ,"CHK %15.8g %15.8g %15.8g %15.8g\n",error0,grad,a,b);
206 //fprintf(ioQQQ,"Pred 1 %15.8g %15.8g %15.8g\n",error1,f1,error0+ft*(grad+ft/(ft-f2)*(b+a*ft)));
207 //fprintf(ioQQQ,"Pred 2 %15.8g %15.8g %15.8g\n",error2,f2,error0+f2*(grad+f2/(ft-f2)*(b+a*f2)));
208 //f1 = (-b+sqrt(b*b-3.*a*grad*(ft-f2)))/(3.*a);
209 //fprintf(ioQQQ,"Pred x %g %g\n",f1,error0+f1*(grad+f1/(ft-f2)*(b+a*f1)));
210 if ( a != 0.0 )
211 {
212 f1 = 1.-3.*(a/b)*(grad/b)*(ft-f2);
213 if (f1 > 0.)
214 f1 = b/(3.*a)*(sqrt(f1)-1.);
215 else
216 f1 = -b/(3.*a);
217 }
218 else
219 {
220 f1 = -grad/(2.*b);
221 }
222 //fprintf(ioQQQ,"Pred n %g %g\n",f1,error0+f1*(grad+f1/(ft-f2)*(b+a*f1)));
223 pred = error0+f1*(grad+f1/(ft-f2)*(b+a*f1));
224 f2 = ft;
225 }
226 error2 = error1;
227 if (f1 > 0.5*f2 || f1 < 0.)
228 f1 = 0.5*f2;
229 else if (f1 < 0.03*f2)
230 f1 = 0.03*f2;
231 }
232
233 /*
234 * b1vec[] is descent direction
235 * new densities in b2vec[]
236 */
237
238 // Count number of attempts required to find new position
239 conv.incrementCounter(NEWTON_LOOP);
240 if (f1 > 1e-6)
241 {
242 for( i=0; i < n; i++ )
243 {
244 b2vec[i] = b0vec[i]-b1vec[i]*f1;
245 }
246 }
247 else
248 {
249 // Try again with larger rlimit
250 lgOK = false;
251 for( i=0; i < n; i++ )
252 {
253 b2vec[i] = b0vec[i];
254 }
255 break;
256 }
257 }
258 if (0 && LOOPMAX == loop)
259 {
260 double rvmax = 0., rval;
261 int imax=0;
262 for( i=0; i < n; ++i )
263 {
264 if (b0vec[i] != 0.)
265 rval = b1vec[i]/b0vec[i];
266 else
267 rval = 0.;
268 fprintf(ioQQQ,"%7s %11.4g %11.4g\n",
269 groupspecies[i]->label.c_str(),rval,b0vec[i]);
270 if (fabs(rval) > rvmax)
271 {
272 rvmax = fabs(rval);
273 imax = i;
274 }
275 }
276 fprintf(ioQQQ,"Biggest is %s\n",groupspecies[imax]->label.c_str());
277 if (0)
278 { // Verify Jacobian
279 long j;
280 for( j=0; j < n; j++ )
281 {
282 b1vec[j] = b0vec[j];
283 }
284 bool lgConserve;
285 jacobn(MoleMap, b1vec,get_ptr(ervals1),get_ptr(amat),false,&lgConserve);
286 for( i=0; i < n; i++ )
287 {
288 for( j=0; j < n; j++ )
289 {
290 b1vec[j] = b0vec[j];
291 }
292 double db = 1e-3*fabs(b0vec[i])+1e-9;
293 b1vec[i] += db;
294 db = b1vec[i]-b0vec[i];
295 jacobn(MoleMap, b1vec,get_ptr(escale),get_ptr(amat),false,&lgConserve);
296 for( j=0; j < n; j++ )
297 {
298 double e1 = MAT(amat,i,j);
299 double e2 = (escale[j]-ervals1[j])/db;
300 if (fabs(e1-e2) > 0.01*fabs(e1+e2))
301 fprintf(ioQQQ,"%7s %7s %11.4g %11.4g %11.4g\n",
302 groupspecies[i]->label.c_str(),groupspecies[j]->label.c_str(),e1,e2,
303 ervals1[j]/db);
304 }
305 }
306 }
307 exit(-1);
308 }
309
310 *error = (realnum) MIN2(error1,1e30);
311 *eqerror = (realnum) MIN2(erroreq,1e30);
312
313 return lgOK;
314}
315
317{
318 if( speciesToPrint==NULL )
319 {
320 fprintf( ioQQQ,"\n NULL species found in mole_print_species_reactions.\n" );
321 return;
322 }
323
324 long numReacts = 0;
325
326 fprintf( ioQQQ,"\n Reactions involving species %s:\n", speciesToPrint->label.c_str() );
327 for( mole_reaction_i p=mole_priv::reactab.begin(); p!=mole_priv::reactab.end(); ++p )
328 {
329 mole_reaction &rate = *p->second;
330 for( long i=0; i<rate.nreactants; i++ )
331 {
332 molecule *sp = rate.reactants[i];
333 if(rate.rvector[i]==NULL && sp==speciesToPrint )
334 {
335 double drate = mole.reaction_rks[ rate.index ];
336 for (long j=0; j<rate.nreactants; ++j)
337 {
338 drate *= mole.species[ rate.reactants[j]->index ].den;
339 }
340 fprintf( ioQQQ,"%s rate = %g\n", rate.label.c_str(), drate );
341 numReacts++;
342 }
343 }
344 for( long i=0; i<rate.nproducts; i++ )
345 {
346 molecule *sp = rate.products[i];
347 if(rate.pvector[i]==NULL && sp==speciesToPrint )
348 {
349 double drate = mole.reaction_rks[ rate.index ];
350 for (long j=0; j<rate.nreactants; ++j)
351 {
352 drate *= mole.species[ rate.reactants[j]->index ].den;
353 }
354 fprintf( ioQQQ,"%s rate = %g\n", rate.label.c_str(), drate );
355 numReacts++;
356 }
357 }
358 }
359 fprintf( ioQQQ," End of reactions involving species %s. There were %li.\n", speciesToPrint->label.c_str(), numReacts );
360
361 return;
362}
363
364STATIC void mole_system_error(long n, long merror,
365 const valarray<double> &a, const valarray<double> &b)
366{
367 fprintf( ioQQQ, " CO_solve getrf_wrapper error %ld",(long int) merror );
368 if( merror > 0 && merror <= n )
369 {
370 fprintf( ioQQQ," - problem with species %s\n\n", groupspecies[merror-1]->label.c_str() );
371 fprintf( ioQQQ,"index \t Row A(i,%li)\t Col A(%li,j) \t B \t Species\n", merror, merror );
372 for( long index=1; index<=n; index++ )
373 {
374 fprintf( ioQQQ,"%li\t%+.4e\t%+.4e\t%+.4e\t%s\n", index,
375 a[(merror-1)*n + index - 1],
376 a[(index -1)*n + merror- 1],
377 b[index-1],
378 groupspecies[index-1]->label.c_str() );
379 }
380
382 }
383
384 fprintf( ioQQQ,"\n" );
385}
386
387int32 solve_system(const valarray<double> &a, valarray<double> &b,
388 long int n, error_print_t error_print)
389{
390 int32 merror;
391 valarray<int32> ipiv(n);
392
393 const int nrefine=3;
394 int i, j, k;
395 valarray<double> lufac(n*n),oldb(n),err(n);
396
397 ASSERT(a.size() == size_t(n*n));
398 ASSERT(b.size() == size_t(n));
399
400 DEBUG_ENTRY("solve_system()");
401
402 lufac = a;
403
404 if (nrefine>0)
405 {
406 oldb = b;
407 }
408
409 merror = 0;
410 getrf_wrapper(n,n,get_ptr(lufac),n,get_ptr(ipiv),&merror);
411 if ( merror != 0 )
412 {
413 if (error_print != NULL)
414 error_print(n,merror,a,b);
415 else
416 fprintf(ioQQQ,"Singular matrix in solve_system\n");
417 return merror; //cdEXIT(EXIT_FAILURE);
418 }
419
420 getrs_wrapper('N',n,1,get_ptr(lufac),n,get_ptr(ipiv),get_ptr(b),n,&merror);
421
422 if ( merror != 0 )
423 {
424 fprintf( ioQQQ, " solve_system: dgetrs finds singular or ill-conditioned matrix\n" );
425 return merror; //cdEXIT(EXIT_FAILURE);
426 }
427
428 for (k=0;k<nrefine;k++)
429 {
430 for (i=0;i<n;i++)
431 {
432 err[i] = oldb[i];
433 }
434 for (j=0;j<n;j++)
435 {
436 for (i=0;i<n;i++)
437 {
438 err[i] -= a[i+j*n]*b[j];
439 }
440 }
441 getrs_wrapper('N',n,1,get_ptr(lufac),n,get_ptr(ipiv),get_ptr(err),n,&merror);
442 if (0)
443 {
444 // Quick-and-dirty condition number estimate
445 // see Golub & Van Loan, 3rd Edn, p128
446 double maxb=0., maxe=0.;
447 for (i=0;i<n;i++)
448 {
449 if (fabs(b[i])>maxb)
450 maxb = fabs(b[i]);
451 if (fabs(err[i])>maxe)
452 maxe = fabs(err[i]);
453 }
454 if (k == 0)
455 {
456 fprintf(ioQQQ,"Max error %g, Condition %g\n",maxe,1e16*maxe/maxb);
457 }
458 }
459 for (i=0;i<n;i++)
460 {
461 b[i] += err[i];
462 }
463 }
464
465 return merror;
466}
static double * amat
FILE * ioQQQ
Definition cddefines.cpp:7
#define ASSERT(exp)
Definition cddefines.h:578
#define STATIC
Definition cddefines.h:97
#define MIN2
Definition cddefines.h:761
float realnum
Definition cddefines.h:103
T * get_ptr(T *v)
Definition cddefines.h:1079
double e2(double t)
Definition service.cpp:299
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
molecule * pvector[MAXPRODUCTS]
Definition mole_priv.h:57
molecule * reactants[MAXREACTANTS]
Definition mole_priv.h:53
molecule * rvector[MAXREACTANTS]
Definition mole_priv.h:54
string label
Definition mole_priv.h:51
molecule * products[MAXPRODUCTS]
Definition mole_priv.h:56
string label
Definition mole.h:142
int index
Definition mole.h:169
t_conv conv
Definition conv.cpp:5
@ NEWTON_LOOP
Definition conv.h:74
@ NEWTON
Definition conv.h:73
t_mole_local mole
Definition mole.cpp:7
#define SMALLABUND
Definition mole.h:13
molecule ** groupspecies
map< string, count_ptr< mole_reaction > >::iterator mole_reaction_i
Definition mole_priv.h:37
@ PRINTSOL
map< string, count_ptr< mole_reaction > > reactab
bool newton_step(GroupMap &MoleMap, const valarray< double > &b0vec, valarray< double > &b2vec, realnum *eqerror, realnum *error, const long n, double *rlimit, double *rmax, valarray< double > &escale, void(*jacobn)(GroupMap &MoleMap, const valarray< double > &b2vec, double *const ervals, double *const amat, const bool lgJac, bool *lgConserve))
int32 solve_system(const valarray< double > &a, valarray< double > &b, long int n, error_print_t error_print)
void mole_print_species_reactions(molecule *speciesToPrint)
#define MAT(a, I_, J_)
STATIC void mole_system_error(long n, long merror, const valarray< double > &a, const valarray< double > &b)
void(* error_print_t)(long, long, const valarray< double > &, const valarray< double > &)
Definition newton_step.h:23
void getrs_wrapper(char trans, long N, long nrhs, double *A, long lda, int32 *ipiv, double *B, long ldb, int32 *info)
void getrf_wrapper(long M, long N, double *A, long lda, int32 *ipiv, int32 *info)