00001
00002
00003
00004
00005
00006
00007
00008 #include "cddefines.h"
00009 #include "newton_step.h"
00010 #include "conv.h"
00011 #include "thirdparty.h"
00012 #include "mole.h"
00013 #include "mole_priv.h"
00014
00015 #define ABSLIM 1e-12
00016 #define ERRLIM 1e-12
00017 # ifdef MAT
00018 # undef MAT
00019 # endif
00020 # define MAT(a,I_,J_) ((a)[(I_)*(n)+(J_)])
00021
00022 STATIC void mole_system_error(long n, long merror,
00023 const valarray<double> &a,
00024 const valarray<double> &b);
00025
00026 enum {PRINTSOL = false};
00027
00028
00029
00030 bool newton_step(GroupMap &MoleMap, const valarray<double> &b0vec, valarray<double> &b2vec,
00031 realnum *eqerror, realnum *error, const long n, double *rlimit, double *rmax,
00032 valarray<double> &escale,
00033 void (*jacobn)(GroupMap &MoleMap,
00034 const valarray<double> &b2vec,
00035 double * const ervals, double * const amat,
00036 const bool lgJac, bool *lgConserve))
00037 {
00038 bool lgOK=true;
00039 long int i,
00040 loop;
00041
00042 double
00043 f1,
00044 f2,
00045 error2,
00046 error1,
00047 error0,
00048 erroreq,
00049 grad,
00050 pred;
00051
00052 valarray<double>
00053 amat(n*n),
00054 b1vec(n),
00055 ervals(n),
00056 ervals1(n);
00057 int32 merror;
00058
00059 DEBUG_ENTRY( "newton_step()" );
00060
00061 for( i=0; i < n; i++ )
00062 {
00063 b2vec[i] = b0vec[i];
00064 }
00065
00066 pred = error0 = error2 = erroreq = grad = f1 = f2 = 0.;
00067
00068 conv.incrementCounter(NEWTON);
00069 const int LOOPMAX = 40;
00070 for (loop=0;loop<LOOPMAX;loop++)
00071 {
00072 bool lgConserve;
00073 jacobn(MoleMap, b2vec,get_ptr(ervals),get_ptr(amat),loop==0,&lgConserve);
00074
00075 if (loop == 0)
00076 {
00077 for( i=0; i < n; i++ )
00078 {
00079 escale[i] = MAT(amat,i,i);
00080 }
00081
00082
00083 *rmax = 0.0;
00084 for( i=0; i<n; ++i )
00085 {
00086 if (-escale[i]>*rmax)
00087 {
00088 *rmax = -escale[i];
00089 }
00090 }
00091 if (*rlimit < 0.0)
00092 {
00093
00094 *rlimit = 1e-19 * (*rmax);
00095 }
00096 else if (*rlimit > *rmax)
00097 {
00098
00099 *rlimit = *rmax;
00100 }
00101 for( i=0; i < n; i++ )
00102 {
00103 if (! lgConserve || ( (! groupspecies[i]->isMonatomic()) || groupspecies[i]->charge < 0 ) )
00104 {
00105
00106
00107 MAT(amat,i,i) -= *rlimit;
00108 }
00109 }
00110 }
00111
00112
00113
00114 erroreq = 0.;
00115 for( i=0; i < n; i++ )
00116 {
00117 double etmp = ervals[i]/(SMALLABUND*(*rmax)+fabs(b0vec[i]*escale[i]));
00118 erroreq += etmp*etmp;
00119 }
00120
00121
00122 for (i=0; i < n; i++)
00123 {
00124 ervals1[i] = ervals[i];
00125 if (! lgConserve || ( (! groupspecies[i]->isMonatomic()) || groupspecies[i]->charge < 0 ) )
00126 ervals1[i] -= (*rlimit)*(b2vec[i]-b0vec[i]);
00127 }
00128
00129 error1 = 0.;
00130 for( i=0; i < n; i++ )
00131 {
00132
00133 double etmp = ervals1[i]/(SMALLABUND*(*rmax)+fabs(b0vec[i]*escale[i]));
00134 etmp *= etmp;
00135 error1 += etmp;
00136 }
00137
00138 if (loop == 0)
00139 {
00140 for( i=0; i < n; i++ )
00141 {
00142 b1vec[i] = ervals1[i];
00143 }
00144 merror = solve_system(amat,b1vec,n,mole_system_error);
00145
00146 if (merror != 0) {
00147 *eqerror = *error = 1e30f;
00148 lgOK = false;
00149 return lgOK;
00150 }
00151 error0 = error1;
00152 grad = -2*error0;
00153 f1 = 1.0;
00154
00155 } else {
00156
00157
00158 if (error1 < (1-2e-4*f1)*error0 || error1 < 1e-20)
00159 {
00160 break;
00161 }
00162
00163 if (loop == 1)
00164 {
00165 f2 = f1;
00166 f1 *= -0.5*f1*grad/(error1-error0-f1*grad);
00167 pred = error0+0.5*grad*f1;
00168 }
00169 else
00170 {
00171 if (0)
00172 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",
00173 loop,f1,error1,erroreq,*rlimit,grad,(error1-error0)/f1,error0,error1-error0,pred);
00174
00175
00176 double a = (error1-error0)/(f1*f1) - (error2-error0)/(f2*f2);
00177 a += grad * (1./f2 - 1./f1);
00178
00179 double b = -f2*(error1-error0)/(f1*f1) + f1*(error2-error0)/(f2*f2);
00180 b += grad * (f2/f1 - f1/f2);
00181 double ft = f1;
00182
00183
00184
00185
00186
00187 if ( a != 0.0 )
00188 {
00189 f1 = 1.-3.*(a/b)*(grad/b)*(ft-f2);
00190 if (f1 > 0.)
00191 f1 = b/(3.*a)*(sqrt(f1)-1.);
00192 else
00193 f1 = -b/(3.*a);
00194 }
00195 else
00196 {
00197 f1 = -grad/(2.*b);
00198 }
00199
00200 pred = error0+f1*(grad+f1/(ft-f2)*(b+a*f1));
00201 f2 = ft;
00202 }
00203 error2 = error1;
00204 if (f1 > 0.5*f2 || f1 < 0.)
00205 f1 = 0.5*f2;
00206 else if (f1 < 0.03*f2)
00207 f1 = 0.03*f2;
00208 }
00209
00210
00211
00212
00213
00214
00215
00216 conv.incrementCounter(NEWTON_LOOP);
00217 if (f1 > 1e-6)
00218 {
00219 for( i=0; i < n; i++ )
00220 {
00221 b2vec[i] = b0vec[i]-b1vec[i]*f1;
00222 }
00223 }
00224 else
00225 {
00226
00227 lgOK = false;
00228 for( i=0; i < n; i++ )
00229 {
00230 b2vec[i] = b0vec[i];
00231 }
00232 break;
00233 }
00234 }
00235 if (0 && LOOPMAX == loop)
00236 {
00237 double rvmax = 0., rval;
00238 int imax=0;
00239 for( i=0; i < n; ++i )
00240 {
00241 if (b0vec[i] != 0.)
00242 rval = b1vec[i]/b0vec[i];
00243 else
00244 rval = 0.;
00245 fprintf(ioQQQ,"%7s %11.4g %11.4g\n",
00246 groupspecies[i]->label.c_str(),rval,b0vec[i]);
00247 if (fabs(rval) > rvmax)
00248 {
00249 rvmax = fabs(rval);
00250 imax = i;
00251 }
00252 }
00253 fprintf(ioQQQ,"Biggest is %s\n",groupspecies[imax]->label.c_str());
00254 if (0)
00255 {
00256 long j;
00257 for( j=0; j < n; j++ )
00258 {
00259 b1vec[j] = b0vec[j];
00260 }
00261 bool lgConserve;
00262 jacobn(MoleMap, b1vec,get_ptr(ervals1),get_ptr(amat),false,&lgConserve);
00263 for( i=0; i < n; i++ )
00264 {
00265 for( j=0; j < n; j++ )
00266 {
00267 b1vec[j] = b0vec[j];
00268 }
00269 double db = 1e-3*fabs(b0vec[i])+1e-9;
00270 b1vec[i] += db;
00271 db = b1vec[i]-b0vec[i];
00272 jacobn(MoleMap, b1vec,get_ptr(escale),get_ptr(amat),false,&lgConserve);
00273 for( j=0; j < n; j++ )
00274 {
00275 double e1 = MAT(amat,i,j);
00276 double e2 = (escale[j]-ervals1[j])/db;
00277 if (fabs(e1-e2) > 0.01*fabs(e1+e2))
00278 fprintf(ioQQQ,"%7s %7s %11.4g %11.4g %11.4g\n",
00279 groupspecies[i]->label.c_str(),groupspecies[j]->label.c_str(),e1,e2,
00280 ervals1[j]/db);
00281 }
00282 }
00283 }
00284 exit(-1);
00285 }
00286
00287 *error = (realnum) MIN2(error1,1e30);
00288 *eqerror = (realnum) MIN2(erroreq,1e30);
00289
00290 return lgOK;
00291 }
00292
00293 void mole_print_species_reactions( molecule *speciesToPrint )
00294 {
00295 if( speciesToPrint==NULL )
00296 {
00297 fprintf( ioQQQ,"\n NULL species found in mole_print_species_reactions.\n" );
00298 return;
00299 }
00300
00301 long numReacts = 0;
00302
00303 fprintf( ioQQQ,"\n Reactions involving species %s:\n", speciesToPrint->label.c_str() );
00304 for( mole_reaction_i p=mole_priv::reactab.begin(); p!=mole_priv::reactab.end(); ++p )
00305 {
00306 mole_reaction &rate = *p->second;
00307 for( long i=0; i<rate.nreactants; i++ )
00308 {
00309 molecule *sp = rate.reactants[i];
00310 if(rate.rvector[i]==NULL && sp==speciesToPrint )
00311 {
00312 double drate = mole.reaction_rks[ rate.index ];
00313 for (long j=0; j<rate.nreactants; ++j)
00314 {
00315 drate *= mole.species[ rate.reactants[j]->index ].den;
00316 }
00317 fprintf( ioQQQ,"%s rate = %g\n", rate.label.c_str(), drate );
00318 numReacts++;
00319 }
00320 }
00321 for( long i=0; i<rate.nproducts; i++ )
00322 {
00323 molecule *sp = rate.products[i];
00324 if(rate.pvector[i]==NULL && sp==speciesToPrint )
00325 {
00326 double drate = mole.reaction_rks[ rate.index ];
00327 for (long j=0; j<rate.nreactants; ++j)
00328 {
00329 drate *= mole.species[ rate.reactants[j]->index ].den;
00330 }
00331 fprintf( ioQQQ,"%s rate = %g\n", rate.label.c_str(), drate );
00332 numReacts++;
00333 }
00334 }
00335 }
00336 fprintf( ioQQQ," End of reactions involving species %s. There were %li.\n", speciesToPrint->label.c_str(), numReacts );
00337
00338 return;
00339 }
00340
00341 STATIC void mole_system_error(long n, long merror,
00342 const valarray<double> &a, const valarray<double> &b)
00343 {
00344 fprintf( ioQQQ, " CO_solve getrf_wrapper error %ld",(long int) merror );
00345 if( merror > 0 && merror <= n )
00346 {
00347 fprintf( ioQQQ," - problem with species %s\n\n", groupspecies[merror-1]->label.c_str() );
00348 fprintf( ioQQQ,"index \t Row A(i,%li)\t Col A(%li,j) \t B \t Species\n", merror, merror );
00349 for( long index=1; index<=n; index++ )
00350 {
00351 fprintf( ioQQQ,"%li\t%+.4e\t%+.4e\t%+.4e\t%s\n", index,
00352 a[(merror-1)*n + index - 1],
00353 a[(index -1)*n + merror- 1],
00354 b[index-1],
00355 groupspecies[index-1]->label.c_str() );
00356 }
00357
00358 mole_print_species_reactions( groupspecies[merror-1] );
00359 }
00360
00361 fprintf( ioQQQ,"\n" );
00362 }
00363
00364 int32 solve_system(const valarray<double> &a, valarray<double> &b,
00365 long int n, error_print_t error_print)
00366 {
00367 int32 merror;
00368 valarray<int32> ipiv(n);
00369
00370 const int nrefine=3;
00371 int i, j, k;
00372 valarray<double> lufac(n*n),oldb(n),err(n);
00373
00374 ASSERT(a.size() == size_t(n*n));
00375 ASSERT(b.size() == size_t(n));
00376
00377 DEBUG_ENTRY("solve_system()");
00378
00379 lufac = a;
00380
00381 if (nrefine>0)
00382 {
00383 oldb = b;
00384 }
00385
00386 merror = 0;
00387 getrf_wrapper(n,n,get_ptr(lufac),n,get_ptr(ipiv),&merror);
00388 if ( merror != 0 )
00389 {
00390 if (error_print != NULL)
00391 error_print(n,merror,a,b);
00392 else
00393 fprintf(ioQQQ,"Singular matrix in solve_system\n");
00394 return merror;
00395 }
00396
00397 getrs_wrapper('N',n,1,get_ptr(lufac),n,get_ptr(ipiv),get_ptr(b),n,&merror);
00398
00399 if ( merror != 0 )
00400 {
00401 fprintf( ioQQQ, " solve_system: dgetrs finds singular or ill-conditioned matrix\n" );
00402 return merror;
00403 }
00404
00405 for (k=0;k<nrefine;k++)
00406 {
00407 for (i=0;i<n;i++)
00408 {
00409 err[i] = oldb[i];
00410 }
00411 for (j=0;j<n;j++)
00412 {
00413 for (i=0;i<n;i++)
00414 {
00415 err[i] -= a[i+j*n]*b[j];
00416 }
00417 }
00418 getrs_wrapper('N',n,1,get_ptr(lufac),n,get_ptr(ipiv),get_ptr(err),n,&merror);
00419 if (0)
00420 {
00421
00422
00423 double maxb=0., maxe=0.;
00424 for (i=0;i<n;i++)
00425 {
00426 if (fabs(b[i])>maxb)
00427 maxb = fabs(b[i]);
00428 if (fabs(err[i])>maxe)
00429 maxe = fabs(err[i]);
00430 }
00431 if (k == 0)
00432 {
00433 fprintf(ioQQQ,"Max error %g, Condition %g\n",maxe,1e16*maxe/maxb);
00434 }
00435 }
00436 for (i=0;i<n;i++)
00437 {
00438 b[i] += err[i];
00439 }
00440 }
00441
00442 return merror;
00443 }