Changeset 2028 for branches/newmole
- Timestamp:
- 05/06/08 16:35:36 (8 months ago)
- Location:
- branches/newmole/source
- Files:
-
- 8 modified
-
molcol.cpp (modified) (3 diffs)
-
mole.h (modified) (1 diff)
-
mole_drive.cpp (modified) (3 diffs)
-
mole_priv.h (modified) (1 diff)
-
mole_solve.cpp (modified) (6 diffs)
-
newton_step.cpp (modified) (11 diffs)
-
newton_step.h (modified) (1 diff)
-
pressure_change.cpp (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
branches/newmole/source/molcol.cpp
r1780 r2028 19 19 20 20 DEBUG_ENTRY( "molcol()" ); 21 22 if( strcmp(chLabel,"PRIN") == 0 )23 {24 /* total hydrogen column density, all forms */25 fprintf( ioMEAN, "\n Log10 Column density (cm^-2)\n");26 fprintf( ioMEAN, " Htot :");27 fprintf( ioMEAN, "%7.3f",log10(MAX2(SMALLFLOAT,colden.colden[ipCOL_HTOT])));28 fprintf( ioMEAN, " HII :");29 fprintf( ioMEAN, "%7.3f",log10(MAX2(SMALLFLOAT,colden.colden[ipCOL_Hp])));30 fprintf( ioMEAN, " HI :");31 fprintf( ioMEAN, "%7.3f",log10(MAX2(SMALLFLOAT,colden.colden[ipCOL_H0])));32 fprintf( ioMEAN, "\n");33 fprintf( ioMEAN, " H- :");34 fprintf( ioMEAN, "%7.3f",log10(MAX2(SMALLFLOAT,colden.colden[ipCOL_HMIN])));35 fprintf( ioMEAN, " H2g :");36 fprintf( ioMEAN, "%7.3f",log10(MAX2(SMALLFLOAT,colden.colden[ipCOL_H2g])));37 fprintf( ioMEAN, " H2* :");38 fprintf( ioMEAN, "%7.3f",log10(MAX2(SMALLFLOAT,colden.colden[ipCOL_H2s])));39 fprintf( ioMEAN, " H2+ :");40 fprintf( ioMEAN, "%7.3f",log10(MAX2(SMALLFLOAT,colden.colden[ipCOL_H2p])));41 fprintf( ioMEAN, " HeH+ :");42 fprintf( ioMEAN, "%7.3f",log10(MAX2(SMALLFLOAT,colden.colden[ipCOL_HeHp] )));43 fprintf( ioMEAN, "\n");44 fprintf( ioMEAN, " H3+ :");45 fprintf( ioMEAN, "%7.3f",log10(MAX2(SMALLFLOAT,colden.colden[ipCOL_H3p] )));46 fprintf( ioMEAN, "\n");47 }48 21 49 22 /* call large H2 and CO column density routines which will do their jobs */ … … 72 45 else if( strcmp(chLabel,"PRIN") == 0 ) 73 46 { 74 /* print the molecular column densities 75 * want to print all the molecules, not including the atoms/ions 76 * that are part of the co solver. use first to print them all */ 77 /*for( i=0; i < mole.num_comole_calc; i++ )*/ 47 /* print the molecular column densities */ 78 48 int j=0; 79 49 for( i=0; i < mole.num_calc; i++ ) … … 84 54 fprintf( ioMEAN, "%7.3f",log10(MAX2(SMALLFLOAT,mole.list[i]->column ))); 85 55 j++; 56 fixit(); /* Separate elemental chem more clearly -- would need to refactor mole.nElem */ 86 57 if( j == 8 ) 87 58 { -
branches/newmole/source/mole.h
r1742 r2028 158 158 159 159 /* Historical solution data */ 160 realnum den_reinit; /* density in first zone of previous iteration */161 160 realnum column_old; /** total column density in previous iteration */ 162 161 }; -
branches/newmole/source/mole_drive.cpp
r2027 r2028 32 32 static const double MOLETOLER = 0.10; 33 33 34 /*lgMolecSetup use previous iteration's converged solution to initialize if available */35 STATIC void MolecSetup (void);36 34 STATIC void mole_effects(void); 37 35 STATIC void mole_h_rate_diagnostics(void); … … 47 45 DEBUG_ENTRY( "mole_drive()" ); 48 46 49 MolecSetup();50 51 47 mole_update_species_cache(); /* Update densities of species controlled outside the chemical network */ 52 48 … … 69 65 70 66 mole_effects(); 71 72 return;73 }74 75 STATIC void MolecSetup(void)76 {77 long int i;78 /* this will be used to rescale old saved abundances in constant pressure model */79 static realnum hden_save_prev_call;80 const bool DEBUG_MOLECAVER = false;81 82 DEBUG_ENTRY( "MolecSetup()" );83 84 /* mole_Init set this to -2 when code initialized, so negative85 * number shows very first call in this model */86 /* >>chng 05 jul 15, do not init if we have never evaluated CO in this iteration87 * and we are below limit where it should be evaluated */88 if( mole.iteration < 0 || true ) /* ???? */89 {90 91 /* very first attempt at ever obtaining a solution -92 * called one time per model since co.iteration_co set negative93 * when cdInit called */94 95 /* >>chng 05 jun 24, during map phase do not reset molecules to zero96 * since we probably have a better estimate right now */97 98 /* we should have a neutral carbon solution at this point */99 //ASSERT( !dense.lgElmtOn[ipCARBON] || dense.xIonDense[ipCARBON][0]>SMALLFLOAT );100 101 /* set iteration flag */102 mole.iteration = iteration;103 mole.nzone = nzone;104 105 /* for constant pressure models when molecules are reset on second106 * and higher iterations, total density will be different, so107 * must take this into account when abundances are reset */108 hden_save_prev_call = dense.gas_phase[ipHYDROGEN];109 110 if( DEBUG_MOLECAVER )111 fprintf(ioQQQ," DEBUG MolecSetup iteration %li zone %li zeroing since very first call. H2/H=%.2e\n",112 iteration,nzone,hmi.H2_total/dense.gas_phase[ipHYDROGEN]);113 }114 else if( iteration > mole.iteration )115 {116 realnum ratio = dense.gas_phase[ipHYDROGEN] / hden_save_prev_call;117 118 /* this is first call on new iteration, reset119 * to first initial abundances from previous iteration */120 for( i=0; i < mole.num_calc; i++ )121 {122 mole.list[i]->den = mole.list[i]->den_reinit*ratio;123 }124 125 mole.iteration = iteration;126 mole.nzone = nzone;127 128 if( DEBUG_MOLECAVER )129 fprintf(ioQQQ," DEBUG MolecSetup iteration %li zone %li resetting since first call on new iteration. H2/H=%.2e\n",130 iteration,131 nzone,132 hmi.H2_total/dense.gas_phase[ipHYDROGEN]);133 }134 else if( iteration == mole.iteration && nzone==mole.nzone+1 )135 {136 /* this branch, second zone with solution, so save previous137 * zone's solution to reset things in next iteration */138 for( i=0; i < mole.num_calc; i++ )139 {140 mole.list[i]->den_reinit = (realnum) mole.list[i]->den;141 }142 143 mole.nzone = -2;144 hden_save_prev_call = dense.gas_phase[ipHYDROGEN];145 146 if( DEBUG_MOLECAVER )147 fprintf(ioQQQ,"DEBUG MolecSetup iteration %li zone %li second zone on new iteration, saving reset.\n", iteration,nzone);148 }149 67 150 68 return; -
branches/newmole/source/mole_priv.h
r1770 r2028 17 17 public: 18 18 double fion[LIMELM][N_MOLE_ION]; 19 void updateMolecules( double *b2);19 void updateMolecules(const double * const b2); 20 20 void setup(double *b0vec); 21 21 }; -
branches/newmole/source/mole_solve.cpp
r1835 r2028 35 35 void check_co_ion_converge(void); 36 36 37 STATIC void funjac(const valarray<double> &b2vec, double * const ervals, double * const amat, const bool lgJac); 38 37 39 STATIC void mole_h_fixup(void); 40 41 #define SMALLABUND 1e-24 38 42 39 43 double mole_solve() … … 42 46 realnum 43 47 error; 44 valarray<double> b2vec(mole.num_compacted); 48 valarray<double> oldmols(mole.num_compacted), 49 newmols(mole.num_compacted); 45 50 46 51 DEBUG_ENTRY( "mole_solve()" ); … … 89 94 nPrevBad = nBad; 90 95 91 MoleMap.setup(&b2vec[0]); 92 newton_step(b2vec, &nBad, &error, mole.num_compacted); 93 MoleMap.updateMolecules( &b2vec[0] ); 96 MoleMap.setup(&oldmols[0]); 97 newton_step(oldmols, newmols, &error, mole.num_compacted, funjac); 98 99 /* check for negative populations */ 100 nBad = 0; 101 long fracneg = 0.; 102 long iworst = -1; 103 for( long i=0; i < mole.num_compacted; i++ ) 104 { 105 if( newmols[i] < 0. ) 106 { 107 if(-newmols[i] > fracneg*SDIV(oldmols[i]) && oldmols[i] > SMALLABUND) 108 { 109 fracneg = -newmols[i]/oldmols[i]; 110 iworst = i; 111 } 112 if (oldmols[i] > SMALLABUND || newmols[i] < -SMALLABUND) 113 ++nBad; 114 } 115 } 116 117 MoleMap.updateMolecules( &newmols[0] ); 94 118 95 119 //fprintf(stdout,"Mole zone %ld -- %7.2f loop %d error %15.8g nBad %d\n",nzone,fnzone,i,error,nBad); … … 248 272 #define ABSLIM 1e-12 249 273 #define ERRLIM 1e-12 250 #define SMALLABUND 1e-24251 274 # ifdef MAT 252 275 # undef MAT … … 260 283 GroupMap MoleMap; 261 284 262 void funjac(valarray<double> &b2vec, double *ervals, double *amat,bool lgJac)285 STATIC void funjac(const valarray<double> &b2vec, double * const ervals, double * const amat, const bool lgJac) 263 286 { 264 287 long … … 520 543 } 521 544 522 void GroupMap::updateMolecules( double *b2)545 void GroupMap::updateMolecules(const double * const b2) 523 546 { 524 547 long int -
branches/newmole/source/newton_step.cpp
r1851 r2028 11 11 #define ABSLIM 1e-12 12 12 #define ERRLIM 1e-12 13 #define SMALLABUND 1e-2414 13 # ifdef MAT 15 14 # undef MAT … … 21 20 enum {PRINTSOL = false}; 22 21 23 extern void funjac(valarray<double> &b2vec, double *ervals, double *amat, bool lgJac);24 25 22 /* mole_newton_step -- improve balance in chemical network along 26 23 * descent direction, step limited to ensure improvement */ 27 void newton_step(valarray<double> &b2vec, int *nBad, realnum *error, long n) 24 void newton_step(const valarray<double> &b0vec, valarray<double> &b2vec, realnum *error, const long n, 25 void (*jacobn)(const valarray<double> &b2vec, double * const ervals, double * const amat, const bool lgJac)) 28 26 { 29 27 long int i, … … 31 29 iworst; 32 30 33 int printsol = PRINTSOL;34 35 31 double 36 32 etmp, 37 33 emax, 38 fracneg,39 34 fstep, 40 35 flast, … … 45 40 valarray<double> 46 41 amat(n*n), 47 b0vec(n),48 42 b1vec(n), 49 43 escale(n), … … 51 45 int32 merror; 52 46 53 DEBUG_ENTRY( " mole_newton_step()" );47 DEBUG_ENTRY( "newton_step()" ); 54 48 55 49 for( i=0; i < n; i++ ) 56 50 { 57 b 0vec[i] = b2vec[i];51 b2vec[i] = b0vec[i]; 58 52 } 59 53 … … 61 55 error1 = grad = 0.; 62 56 fstep = flast = 1.; 63 *nBad = 0;64 57 65 58 for (loop=0;loop<20;loop++) … … 67 60 flast = fstep; 68 61 69 funjac(b2vec,&ervals[0],&amat[0],loop==0);62 jacobn(b2vec,&ervals[0],&amat[0],loop==0); 70 63 71 64 if (loop == 0) … … 89 82 iworst = i; 90 83 } 91 92 if(printsol)93 fprintf(ioQQQ,"%15.7g\t",etmp);94 84 error0 += etmp; 95 85 } 96 97 // fprintf(stdout,"Loop %ld error %15.8g fstep %15.8g; worst species %s abund %g error %g\n",98 // loop,error0,fstep,groupspecies[iworst]->label,groupspecies[iworst]->den,,emax);99 86 100 87 if (loop == 0) … … 107 94 error1 = error0; 108 95 grad = -2*error1; 109 /* b1vec is now descent direction */110 111 //fprintf(ioQQQ,"Step factor %g bad species %s abund %g step %g\n",112 // fstep,iworst!=-1?groupspecies[iworst]->label:"None",113 // iworst!=-1?b0vec[iworst]:0.,114 // iworst!=-1?-b1vec[iworst]:0.);115 96 } else { 116 //fprintf(ioQQQ,"Step factor %15.8g bad species %d %s abund %15.8g step %15.8g\n",117 // fstep,iworst,iworst!=-1?groupspecies[iworst]->label:"None",118 // iworst!=-1?b0vec[iworst]:0.,119 // iworst!=-1?-b1vec[iworst]:0.);120 //fprintf(stdout,"Grads %15.8g %15.8g (%15.8g %15.8g)\n",grad,(error0-error1)/fstep,error0,error1);121 97 if (error0 < (1-2e-4*fstep)*error1 || error0 < 1e-20) 122 98 { … … 124 100 } 125 101 fstep = -0.5*grad*fstep*fstep/SDIV(error0-error1-grad*fstep); 126 //fprintf(stdout,"Errs1 loop %ld %g now %g fac %g was %g pred %g\n",loop,127 // error1,error0,fstep,flast,128 // error1+grad*fstep+(fstep*fstep/(flast*flast))*(error0-(error1+grad*flast)));129 102 if (fstep > 0.9*flast || fstep < 0.) 130 103 fstep = 0.5*flast; … … 133 106 } 134 107 135 //if (iworst != -1)136 // fprintf(stdout,"Loop %ld, worst error for %s: %g of %g, step %g of %g %g\n",137 // loop, groupspecies[iworst]->label,emax,error0,-b1vec[iworst]*fstep,b0vec[iworst],fstep);138 139 //fprintf(ioQQQ,"Loops %ld %g %g\n",loop,error0,(1-2e-4*fstep)*error1);140 141 108 /* 142 109 * b1vec[] is descent direction 143 * new densit yin b2vec[]110 * new densities in b2vec[] 144 111 */ 145 112 146 113 for( i=0; i < n; i++ ) 147 114 { 148 115 b2vec[i] = b0vec[i]-b1vec[i]*fstep; 149 116 } 150 151 /* check for negative populations */152 *nBad = 0;153 fracneg = 0.;154 iworst = -1;155 for( i=0; i < n; i++ )156 {157 if( b2vec[i] < 0. )158 {159 // fprintf(stdout,"Species %s negative at %g from %g\n",groupspecies[i]->label,b2vec[i],b0vec[i]);160 /* largest relative change in solution for neg soln */161 /* this can only occur for negative solutions since fracneg starts162 * as zero */163 if(-b2vec[i] > fracneg*SDIV(b0vec[i]) && b0vec[i] > SMALLABUND)164 {165 fracneg = -b2vec[i]/SDIV(b0vec[i]);166 iworst = i;167 }168 if (b0vec[i] > SMALLABUND || b2vec[i] < -SMALLABUND)169 ++*nBad;170 }171 }172 173 /* Expand species groups into mole.list[]->den -- either for174 * re-evaluation of network error, or for communication with175 * outside code */176 177 117 } 178 //fprintf(ioQQQ,"Final fstep %g\n",fstep);179 118 180 119 *error = (realnum) MIN2(error0,1e30); -
branches/newmole/source/newton_step.h
r1770 r2028 10 10 \param *error 11 11 */ 12 void newton_step(valarray<double> &b2vec, int *nBad, realnum *error, long n); 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)); 13 14 14 15 #endif /* _NEWTON_STEP_H_ */ -
branches/newmole/source/pressure_change.cpp
r1942 r2028 274 274 TempChange(phycon.te , false); 275 275 276 /* molecules done in hmole, only change pressure, not abundances */277 276 hmi.H2_total *= (realnum)pressure_change_factor; 278 277 h2.ortho_density *= (realnum)pressure_change_factor; 279 278 h2.para_density *= (realnum)pressure_change_factor; 280 /* molecules done in comole */281 /* >>chng 03 sep 15, upper limit had not included the C, O atoms/ions */282 /*for( ion=0; ion < NUM_HEAVY_MOLEC; ion++ )*/283 279 for( mol=0; mol < mole.num_total; mol++ ) 284 280 {
