cloudy  trunk
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
newton_step.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2022 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 #include "trace.h"
15 #include "save.h"
16 
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  static ConvergenceCounter cctr=conv.register_("NEWTON");
72  ++cctr;
73  const int LOOPMAX = 40;
74  for (loop=0;loop<LOOPMAX;loop++)
75  {
76  bool lgConserve;
77  jacobn(MoleMap, b2vec,get_ptr(ervals),get_ptr(amat),loop==0,&lgConserve);
78 
79  if (loop == 0)
80  {
81  for( i=0; i < n; i++ )
82  {
83  escale[i] = MAT(amat,i,i);
84  }
85  // First time through this case, set rlimit by guesswork based on
86  // maintaining matrix condition
87  *rmax = 0.0;
88  for( i=0; i<n; ++i )
89  {
90  if (-escale[i]>*rmax)
91  {
92  *rmax = -escale[i];
93  }
94  }
95  if (*rlimit < 0.0)
96  {
97  // large generally more stable, small more quickly convergent
98  *rlimit = 1e-19 * (*rmax);
99  }
100  else if (*rlimit > *rmax)
101  {
102  // large generally more stable, small more quickly convergent
103  *rlimit = *rmax;
104  }
105  for( i=0; i < n; i++ )
106  {
107  if (! lgConserve || ( (! groupspecies[i]->isMonatomic()) || groupspecies[i]->charge < 0 ) )
108  {
109  // Only apply *rlimit to rows which haven't been
110  // overwritten by conservation constraints
111  MAT(amat,i,i) -= *rlimit;
112  }
113  }
114  }
115 
116  // External error limit based on difference to state relaxed to
117  // equilibrium
118  erroreq = 0.;
119  for( i=0; i < n; i++ )
120  {
121  double etmp = ervals[i]/(SMALLABUND*(*rmax)+fabs(b0vec[i]*escale[i]));
122  erroreq += etmp*etmp;
123  }
124 
125  if (trace.lgTrace)
126  {
127  fprintf(ioQQQ,"Newton step loop %ld error %15.8g\n",loop,erroreq);
128  }
129 
130  // Internal (step) limit based on finite-(pseudo-)timestep solution
131  for (i=0; i < n; i++)
132  {
133  ervals1[i] = ervals[i];
134  if (! lgConserve || ( (! groupspecies[i]->isMonatomic()) || groupspecies[i]->charge < 0 ) )
135  ervals1[i] -= (*rlimit)*(b2vec[i]-b0vec[i]);
136  }
137 
138  error1 = 0.;
139  long maxi = -1;
140  double emax0 = 0.0, emax1 = 0.0;
141  for( i=0; i < n; i++ )
142  {
143  if (! lgConserve || ( (! groupspecies[i]->isMonatomic()) || groupspecies[i]->charge < 0 ) )
144  {
145  /* Scale the errors weighting to ensure trace species are accurate */
146  double etmp = ervals1[i]/(SMALLABUND*(*rmax)+fabs(b0vec[i]*escale[i]));
147  double etmp0 = ervals0[i]/(SMALLABUND*(*rmax)+fabs(b0vec[i]*escale[i]));
148  etmp *= etmp;
149  etmp0 *= etmp0;
150  if (fabs(etmp-etmp0) > fabs(emax1-emax0) )
151  {
152  maxi = i;
153  emax1 = etmp;
154  emax0 = etmp0;
155  }
156  error1 += etmp;
157  }
158  }
159 
160  if (loop == 0)
161  {
162  for( i=0; i < n; i++ )
163  {
164  ervals0[i] = ervals1[i];
165  b1vec[i] = ervals1[i];
166  }
167  merror = solve_system(amat,b1vec,n,mole_system_error);
168 
169  if (merror != 0) {
170  *eqerror = *error = 1e30f;
171  lgOK = false;
172  return lgOK;
173  }
174  error0 = error1;
175  grad = -2*error0;
176  f1 = 1.0;
177 
178  } else {
179  //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",
180  // loop,f1,error1,erroreq,*rlimit,grad,(error1-error0)/f1,error0,error1-error0);
181  if (error1 < (1-2e-4*f1)*error0 || error1 < 1e-20)
182  {
183  break;
184  }
185  // Backtrack using quadratic or cubic fit, ref Dennis & Schnabel 1996
186  if (loop == 1)
187  {
188  f2 = f1;
189  f1 *= -0.5*f1*grad/(error1-error0-f1*grad);
190  pred = error0+0.5*grad*f1;
191  }
192  else
193  {
194  if (0)
195  {
196  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",
197  loop,f1,error1,erroreq,*rlimit,grad,(error1-error0)/f1,error0,error1-error0,pred
198  );
199  if (maxi != -1)
200  fprintf(ioQQQ,"Maxi %ld %s emax from %11.4g of %11.4g -> %11.4g of %11.4g, diff %11.4g\n",
201  maxi,groupspecies[maxi]->label.c_str(),emax0,error0, emax1,error1,emax1-emax0);
202  }
203  // NB we must be careful how a is calculated because terms can be nearly equal and/or
204  // vastly different from other terms.
205  double a = (error1-error0)/(f1*f1) - (error2-error0)/(f2*f2);
206  a += grad * (1./f2 - 1./f1);
207  // similarly calculate b
208  double b = -f2*(error1-error0)/(f1*f1) + f1*(error2-error0)/(f2*f2);
209  b += grad * (f2/f1 - f1/f2);
210  double ft = f1;
211  //fprintf(ioQQQ,"CHK %15.8g %15.8g %15.8g %15.8g\n",error0,grad,a,b);
212  //fprintf(ioQQQ,"Pred 1 %15.8g %15.8g %15.8g\n",error1,f1,error0+ft*(grad+ft/(ft-f2)*(b+a*ft)));
213  //fprintf(ioQQQ,"Pred 2 %15.8g %15.8g %15.8g\n",error2,f2,error0+f2*(grad+f2/(ft-f2)*(b+a*f2)));
214  //f1 = (-b+sqrt(b*b-3.*a*grad*(ft-f2)))/(3.*a);
215  //fprintf(ioQQQ,"Pred x %g %g\n",f1,error0+f1*(grad+f1/(ft-f2)*(b+a*f1)));
216  if ( a != 0.0 )
217  {
218  f1 = 1.-3.*(a/b)*(grad/b)*(ft-f2);
219  if (f1 > 0.)
220  f1 = b/(3.*a)*(sqrt(f1)-1.);
221  else
222  f1 = -b/(3.*a);
223  }
224  else
225  {
226  f1 = -grad/(2.*b);
227  }
228  //fprintf(ioQQQ,"Pred n %g %g\n",f1,error0+f1*(grad+f1/(ft-f2)*(b+a*f1)));
229  pred = error0+f1*(grad+f1/(ft-f2)*(b+a*f1));
230  f2 = ft;
231  }
232  error2 = error1;
233  if (f1 > 0.5*f2 || f1 < 0.)
234  f1 = 0.5*f2;
235  else if (f1 < 0.03*f2)
236  f1 = 0.03*f2;
237  }
238 
239  /*
240  * b1vec[] is descent direction
241  * new densities in b2vec[]
242  */
243 
244  // Count number of attempts required to find new position
245  static ConvergenceCounter cctrl=conv.register_("NEWTON_LOOP");
246  ++cctrl;
247  if (f1 > 1e-6)
248  {
249  for( i=0; i < n; i++ )
250  {
251  b2vec[i] = b0vec[i]-b1vec[i]*f1;
252  }
253  }
254  else
255  {
256  // Try again with larger rlimit
257  lgOK = false;
258  for( i=0; i < n; i++ )
259  {
260  b2vec[i] = b0vec[i];
261  }
262  break;
263  }
264  }
265  if (0 && LOOPMAX == loop)
266  {
267  double rvmax = 0., rval;
268  int imax=0;
269  for( i=0; i < n; ++i )
270  {
271  if (b0vec[i] != 0.)
272  rval = b1vec[i]/b0vec[i];
273  else
274  rval = 0.;
275  fprintf(ioQQQ,"%7s %11.4g %11.4g\n",
276  groupspecies[i]->label.c_str(),rval,b0vec[i]);
277  if (fabs(rval) > rvmax)
278  {
279  rvmax = fabs(rval);
280  imax = i;
281  }
282  }
283  fprintf(ioQQQ,"Biggest is %s\n",groupspecies[imax]->label.c_str());
284  if (0)
285  { // Verify Jacobian
286  long j;
287  for( j=0; j < n; j++ )
288  {
289  b1vec[j] = b0vec[j];
290  }
291  bool lgConserve;
292  jacobn(MoleMap, b1vec,get_ptr(ervals1),get_ptr(amat),false,&lgConserve);
293  for( i=0; i < n; i++ )
294  {
295  for( j=0; j < n; j++ )
296  {
297  b1vec[j] = b0vec[j];
298  }
299  double db = 1e-3*fabs(b0vec[i])+1e-9;
300  b1vec[i] += db;
301  db = b1vec[i]-b0vec[i];
302  jacobn(MoleMap, b1vec,get_ptr(escale),get_ptr(amat),false,&lgConserve);
303  for( j=0; j < n; j++ )
304  {
305  double e1 = MAT(amat,i,j);
306  double e2 = (escale[j]-ervals1[j])/db;
307  if (fabs(e1-e2) > 0.01*fabs(e1+e2))
308  fprintf(ioQQQ,"%7s %7s %11.4g %11.4g %11.4g\n",
309  groupspecies[i]->label.c_str(),groupspecies[j]->label.c_str(),e1,e2,
310  ervals1[j]/db);
311  }
312  }
313  }
314  exit(-1);
315  }
316 
317  *error = (realnum) MIN2(error1,1e30);
318  *eqerror = (realnum) MIN2(erroreq,1e30);
319 
320  return lgOK;
321 }
322 
323 STATIC void mole_system_error(long n, long merror,
324  const valarray<double> &a, const valarray<double> &b)
325 {
326  fprintf( ioQQQ, " CO_solve getrf_wrapper error %ld",(long int) merror );
327  if( merror > 0 && merror <= n )
328  {
329  fprintf( ioQQQ," - problem with species %s\n\n", groupspecies[merror-1]->label.c_str() );
330  fprintf( ioQQQ,"index \t Row A(i,%li)\t Col A(%li,j) \t B \t Species\n", merror, merror );
331  for( long index=1; index<=n; index++ )
332  {
333  fprintf( ioQQQ,"%li\t%+.4e\t%+.4e\t%+.4e\t%s\n", index,
334  a[(merror-1)*n + index - 1],
335  a[(index -1)*n + merror- 1],
336  b[index-1],
337  groupspecies[index-1]->label.c_str() );
338  }
339 
341  }
342 
343  fprintf( ioQQQ,"\n" );
344 }
345 
346 int32 solve_system(const valarray<double> &a, valarray<double> &b,
347  long int n, error_print_t error_print)
348 {
349  int32 merror;
350  valarray<int32> ipiv(n);
351 
352  const int nrefine=3;
353  int i, j, k;
354  valarray<double> lufac(n*n),oldb(n),err(n);
355 
356  ASSERT(a.size() == size_t(n*n));
357  ASSERT(b.size() == size_t(n));
358 
359  DEBUG_ENTRY( "solve_system()" );
360 
361  lufac = a;
362 
363  if (nrefine>0)
364  {
365  oldb = b;
366  }
367 
368  merror = 0;
369  getrf_wrapper(n,n,get_ptr(lufac),n,get_ptr(ipiv),&merror);
370  if ( merror != 0 )
371  {
372  if (error_print != NULL)
373  error_print(n,merror,a,b);
374  else
375  fprintf(ioQQQ," PROBLEM singular matrix in solve_system\n");
376  return merror; //cdEXIT(EXIT_FAILURE);
377  }
378 
379  getrs_wrapper('N',n,1,get_ptr(lufac),n,get_ptr(ipiv),get_ptr(b),n,&merror);
380 
381  if ( merror != 0 )
382  {
383  fprintf( ioQQQ, " solve_system: dgetrs finds singular or ill-conditioned matrix\n" );
384  return merror; //cdEXIT(EXIT_FAILURE);
385  }
386 
387  for (k=0;k<nrefine;k++)
388  {
389  for (i=0;i<n;i++)
390  {
391  err[i] = oldb[i];
392  }
393  for (j=0;j<n;j++)
394  {
395  for (i=0;i<n;i++)
396  {
397  err[i] -= a[i+j*n]*b[j];
398  }
399  }
400  getrs_wrapper('N',n,1,get_ptr(lufac),n,get_ptr(ipiv),get_ptr(err),n,&merror);
401  if (0)
402  {
403  // Quick-and-dirty condition number estimate
404  // see Golub & Van Loan, 3rd Edn, p128
405  double maxb=0., maxe=0.;
406  for (i=0;i<n;i++)
407  {
408  if (fabs(b[i])>maxb)
409  maxb = fabs(b[i]);
410  if (fabs(err[i])>maxe)
411  maxe = fabs(err[i]);
412  }
413  if (k == 0)
414  {
415  fprintf(ioQQQ,"Max error %g, Condition %g\n",maxe,1e16*maxe/maxb);
416  }
417  }
418  for (i=0;i<n;i++)
419  {
420  b[i] += err[i];
421  }
422  }
423 
424  return merror;
425 }
void getrs_wrapper(char trans, long N, long nrhs, double *A, long lda, int32 *ipiv, double *B, long ldb, int32 *info)
T * get_ptr(T *v)
Definition: cddefines.h:1109
double e1(double x)
ConvergenceCounter register_(const string name)
Definition: conv.cpp:87
t_conv conv
Definition: conv.cpp:5
#define SMALLABUND
Definition: mole.h:17
FILE * ioQQQ
Definition: cddefines.cpp:7
void(* error_print_t)(long, long, const valarray< double > &, const valarray< double > &)
Definition: newton_step.h:22
#define MIN2(a, b)
Definition: cddefines.h:803
t_trace trace
Definition: trace.cpp:5
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
void getrf_wrapper(long M, long N, double *A, long lda, int32 *ipiv, int32 *info)
#define STATIC
Definition: cddefines.h:118
bool lgTrace
Definition: trace.h:12
float realnum
Definition: cddefines.h:124
int32 solve_system(const valarray< double > &a, valarray< double > &b, long int n, error_print_t error_print)
#define ASSERT(exp)
Definition: cddefines.h:613
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:723
int fprintf(const Output &stream, const char *format,...)
Definition: service.cpp:1121
vector< molecule * > groupspecies
double e2(double x)
#define MAT(a, I_, J_)
Definition: newton_step.cpp:20
STATIC void mole_system_error(long n, long merror, const valarray< double > &a, const valarray< double > &b)
void mole_print_species_reactions(molecule *speciesToPrint)