Changeset 2123 for branches/newmole
- Timestamp:
- 06/24/08 16:01:19 (7 months ago)
- Location:
- branches/newmole/source
- Files:
-
- 3 modified
-
mole_solve.cpp (modified) (7 diffs)
-
newton_step.cpp (modified) (5 diffs)
-
newton_step.h (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
branches/newmole/source/mole_solve.cpp
r2118 r2123 43 43 44 44 #define SMALLABUND 1e-24 45 #ifdef MINPACK46 #include <minpack.h>47 STATIC void fcn(const int *n, const double *x, double *fvec, double *fjac,48 const int *ldfjac, int *iflag );49 #endif50 45 51 46 static double molnum[LIMELM]; … … 73 68 74 69 /* will test against this error for quitting */ 75 # define BIGERROR 1e- 470 # define BIGERROR 1e-8 76 71 /* >>chng 04 jul 20, upper limit had been 6, why? change to 20 77 72 * to get closer to soln */ … … 81 76 nBad = nPrevBad = 0; 82 77 error = -1.; 78 79 valarray<double> escale(mole.num_compacted); 80 double rlimit=-1.; 81 83 82 for(i=0; i < LIM_LOOP;i++) 84 83 { … … 105 104 106 105 MoleMap.setup(&oldmols[0]); 107 if (1)108 {109 newton_step(oldmols, newmols, &error, mole.num_compacted, funjac);110 }111 else112 {113 #ifdef MINPACK114 int info, lwa = (mole.num_compacted*(mole.num_compacted+13))/2;115 valarray<double> fvec(mole.num_compacted),116 fjac(mole.num_compacted*mole.num_compacted),117 wa(lwa);118 double tol = 1e-9;119 for ( long i=0; i < mole.num_compacted; i++ )120 {121 newmols[i] = oldmols[i];122 }123 hybrj1_(fcn, &mole.num_compacted, &newmols[0], &fvec[0], &fjac[0],124 &mole.num_compacted, &tol, &info, &wa[0], &lwa);125 //fprintf(ioQQQ,"hybrj1_ returns info %d\n",info);126 #else127 ;128 #endif129 }130 131 /* check for negative populations */132 nBad = 0;133 double fracneg = 0.;134 long iworst = -1;135 for( long i=0; i < mole.num_compacted; i++ )136 {137 if( newmols[i] < 0. )138 {139 if(-newmols[i] > fracneg*SDIV(oldmols[i]) && oldmols[i] > SMALLABUND)140 {141 fracneg = -newmols[i]/oldmols[i];142 iworst = i;143 }144 if (oldmols[i] > SMALLABUND || newmols[i] < -SMALLABUND)145 ++nBad;146 }147 }148 106 107 long j; 108 109 for (j=0; j<10; j++) 110 { 111 newton_step(oldmols, newmols, &error, mole.num_compacted, &rlimit, escale, funjac); 112 113 /* check for negative populations */ 114 nBad = 0; 115 double fracneg = 0.; 116 long iworst = -1; 117 for( long mol=0; mol < mole.num_compacted; mol++ ) 118 { 119 if( newmols[mol] < 0. ) 120 { 121 if(-newmols[mol] > fracneg*SDIV(oldmols[mol]) && oldmols[mol] > SMALLABUND) 122 { 123 fracneg = -newmols[mol]/oldmols[mol]; 124 iworst = mol; 125 } 126 if (oldmols[mol] > SMALLABUND || newmols[mol] < -SMALLABUND) 127 ++nBad; 128 } 129 } 130 if (0 == nBad) 131 { 132 //fprintf(ioQQQ,"OK at z %ld l %d j %ld t %g e %g\n",nzone,i,j,1./rlimit,error); 133 break; 134 } 135 rlimit = 10.*rlimit; 136 } 137 138 // If it got through with no problems, extend the time limit 139 if (0 == j) 140 rlimit = 0.5*rlimit; 141 149 142 MoleMap.updateMolecules( newmols ); 150 143 … … 172 165 ConvFail( "pops" , "Mole "); 173 166 } 167 168 //fprintf(stdout,"Mole zone %ld -- %7.2f error %15.8g nBad %d loop %d\n",nzone,fnzone,error,nBad,i); 174 169 175 170 mole_return_cached_species(); … … 315 310 GroupMap MoleMap; 316 311 317 #ifdef MINPACK318 STATIC void fcn(const int *n, const double *x, double *fvec, double *fjac,319 const int *ldfjac, int *iflag )320 {321 /* Function to interface funjac with MINPACK routines */322 ASSERT (*ldfjac == *n);323 valarray<double> b2vec(*n), tmpvec(*n);324 for (long i=0; i < *n; i++)325 {326 b2vec[i] = x[i];327 }328 if (*iflag == 1)329 {330 funjac(b2vec, fvec, &tmpvec[0], false);331 }332 else if (*iflag == 2)333 {334 funjac(b2vec, &tmpvec[0], &fjac[0], true);335 }336 }337 #endif338 339 312 STATIC void funjac(const valarray<double> &b2vec, double * const ervals, double * const amat, const bool lgJac) 340 313 { … … 363 336 // Without pseudo-timestepping, need to apply conservation 364 337 // constraint to maintain non-singular matrix 365 lgConserve = false;338 lgConserve = true; 366 339 } 367 340 -
branches/newmole/source/newton_step.cpp
r2118 r2123 25 25 * descent direction, step limited to ensure improvement */ 26 26 void newton_step(const valarray<double> &b0vec, valarray<double> &b2vec, 27 realnum *error, const long n, 27 realnum *error, const long n, double *rlimit, valarray<double> &escale, 28 28 void (*jacobn)(const valarray<double> &b2vec, 29 29 double * const ervals, double * const amat, … … 45 45 erroreq, 46 46 grad, 47 pred, 48 errstp; 47 pred; 49 48 50 49 valarray<double> 51 50 amat(n*n), 52 51 b1vec(n), 53 escale(n),54 52 ervals(n), 55 53 ervals1(n); … … 64 62 65 63 pred = error0 = error2 = erroreq = grad = f1 = f2 = 0.; 66 errstp = 0.0;67 64 68 65 for (loop=0;loop<40;loop++) … … 72 69 if (loop == 0) 73 70 { 74 // Limit timestep to value where linearizations are reasonable75 71 iworst = -1; 76 72 for( i=0; i < n; i++ ) 77 73 { 78 74 escale[i] = MAT(amat,i,i); 79 if (b0vec[i] > 0 && errstp*b0vec[i] < fabs(ervals[i])) 80 { 81 errstp = fabs(ervals[i])/(1e-24+b0vec[i]); 82 iworst = i; 83 } 84 } 85 errstp *= 3.; 86 for (i=0; i < n; i++) 87 { 88 MAT(amat,i,i) -= errstp; 89 } 90 //if (iworst != -1) 91 // fprintf(stdout,"Timescale %g rate %g worst %s %g %g %g\n", 92 // 1./SDIV(errstp),errstp,groupspecies[iworst]->label, 93 // ervals[iworst],b0vec[iworst],escale[iworst]); 94 } 95 96 // Error from fully relaxed equilibrium solution 75 } 76 // First time through this case, set rlimit by guesswork based on 77 // maintaining matrix condition 78 if (*rlimit < 0.0) 79 { 80 *rlimit=0.0; 81 for( i=0; i<n; ++i ) 82 { 83 if (-escale[i]>*rlimit) 84 { 85 *rlimit = -escale[i]; 86 } 87 } 88 *rlimit *= 1e-11; 89 } 90 for( i=0; i < n; i++ ) 91 { 92 MAT(amat,i,i) -= *rlimit; 93 } 94 } 95 96 // External error limit based on difference to state relaxed to (quasi-)equilibrium 97 97 erroreq = 0.; 98 98 for( i=0; i < n; i++ ) 99 99 { 100 etmpeq = ervals[i] ;//(1e-24+fabs(b0vec[i]*escale[i]));100 etmpeq = ervals[i]/(1e-24+fabs(b0vec[i]*escale[i])); 101 101 erroreq += etmpeq*etmpeq; 102 102 } 103 103 104 // Error from timestep-limitedsolution104 // Internal (step) limit based on finite-(pseudo-)timestep solution 105 105 for (i=0; i < n; i++) 106 106 { 107 ervals1[i] = ervals[i]- errstp*(b2vec[i]-b0vec[i]);107 ervals1[i] = ervals[i]-(*rlimit)*(b2vec[i]-b0vec[i]); 108 108 } 109 109 error1 = emax = 0.; … … 111 111 for( i=0; i < n; i++ ) 112 112 { 113 /* S mooth the error mode tailoff*/114 etmp = ervals1[i] ;//(1e-24+fabs(b0vec[i]*escale[i]));113 /* Scale the errors weighting to ensure trace species are accurate */ 114 etmp = ervals1[i]/(1e-24+fabs(b0vec[i]*escale[i])); 115 115 etmp *= etmp; 116 116 if (etmp > emax) -
branches/newmole/source/newton_step.h
r2028 r2123 10 10 \param *error 11 11 */ 12 void newton_step(const valarray<double> &oldmols, valarray<double> &newmols, realnum *error, const long n, 13 void (*jacobn)(const valarray<double> &b2vec, double * const ervals, double * const amat, const bool lgJac)); 12 void newton_step(const valarray<double> &oldmols, valarray<double> &newmols, realnum *error, const long n, 13 double *rlimit, valarray<double> &escale, 14 void (*jacobn)(const valarray<double> &b2vec, double * const ervals, double * const amat, const bool lgJac)); 14 15 15 16 #endif /* _NEWTON_STEP_H_ */
