cloudy  trunk
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 
22 STATIC void mole_system_error(long n, long merror,
23  const valarray<double> &a,
24  const valarray<double> &b);
25 
26 enum {PRINTSOL = false};
27 
28 /* mole_newton_step -- improve balance in chemical network along
29  * descent direction, step limited to ensure improvement */
30 bool 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 
316 void mole_print_species_reactions( molecule *speciesToPrint )
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 
364 STATIC 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 
381  mole_print_species_reactions( groupspecies[merror-1] );
382  }
383 
384  fprintf( ioQQQ,"\n" );
385 }
386 
387 int32 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 }
mole_system_error
STATIC void mole_system_error(long n, long merror, const valarray< double > &a, const valarray< double > &b)
Definition: newton_step.cpp:364
mole_reaction::nreactants
int nreactants
Definition: mole_priv.h:52
SMALLABUND
#define SMALLABUND
Definition: mole.h:13
ioQQQ
FILE * ioQQQ
Definition: cddefines.cpp:7
newton_step.h
realnum
float realnum
Definition: cddefines.h:103
conv.h
STATIC
#define STATIC
Definition: cddefines.h:97
mole.h
mole_reaction::pvector
molecule * pvector[MAXPRODUCTS]
Definition: mole_priv.h:57
molecule::index
int index
Definition: mole.h:169
mole_reaction
Definition: mole_priv.h:49
thirdparty.h
MAT
#define MAT(a, I_, J_)
Definition: newton_step.cpp:20
molecule::charge
int charge
Definition: mole.h:144
getrs_wrapper
void getrs_wrapper(char trans, long N, long nrhs, double *A, long lda, int32 *ipiv, double *B, long ldb, int32 *info)
Definition: thirdparty_lapack.cpp:69
ASSERT
#define ASSERT(exp)
Definition: cddefines.h:578
MIN2
#define MIN2
Definition: cddefines.h:761
PRINTSOL
@ PRINTSOL
Definition: newton_step.cpp:26
mole_priv::reactab
map< string, count_ptr< mole_reaction > > reactab
Definition: mole_species.cpp:60
GroupMap
Definition: mole_priv.h:21
mole
t_mole_local mole
Definition: mole.cpp:7
error_print_t
void(* error_print_t)(long, long, const valarray< double > &, const valarray< double > &)
Definition: newton_step.h:23
cddefines.h
mole_reaction::index
long index
Definition: mole_priv.h:62
e2
double e2(double t)
Definition: service.cpp:299
NEWTON
@ NEWTON
Definition: conv.h:73
mole_reaction::nproducts
int nproducts
Definition: mole_priv.h:52
groupspecies
molecule ** groupspecies
Definition: mole_species.cpp:72
getrf_wrapper
void getrf_wrapper(long M, long N, double *A, long lda, int32 *ipiv, int32 *info)
Definition: thirdparty_lapack.cpp:47
mole_print_species_reactions
void mole_print_species_reactions(molecule *speciesToPrint)
Definition: newton_step.cpp:316
mole_reaction::reactants
molecule * reactants[MAXREACTANTS]
Definition: mole_priv.h:53
mole_reaction::products
molecule * products[MAXPRODUCTS]
Definition: mole_priv.h:56
mole_reaction::label
string label
Definition: mole_priv.h:51
NEWTON_LOOP
@ NEWTON_LOOP
Definition: conv.h:74
newton_step
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))
Definition: newton_step.cpp:30
conv
t_conv conv
Definition: conv.cpp:5
get_ptr
T * get_ptr(T *v)
Definition: cddefines.h:1079
t_conv::incrementCounter
void incrementCounter(const counter_type type)
Definition: conv.h:308
molecule
Definition: mole.h:132
t_mole_local::species
valarray< class molezone > species
Definition: mole.h:398
mole_reaction_i
map< string, count_ptr< mole_reaction > >::iterator mole_reaction_i
Definition: mole_priv.h:37
molecule::label
string label
Definition: mole.h:142
t_mole_local::reaction_rks
vector< double > reaction_rks
Definition: mole.h:400
mole_priv.h
amat
static double * amat
Definition: atom_feii.cpp:173
solve_system
int32 solve_system(const valarray< double > &a, valarray< double > &b, long int n, error_print_t error_print)
Definition: newton_step.cpp:387
DEBUG_ENTRY
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:684
mole_reaction::rvector
molecule * rvector[MAXREACTANTS]
Definition: mole_priv.h:54
Generated by   doxygen 1.8.17