| 1 | Index: conv_temp_eden_ioniz.cpp |
|---|
| 2 | =================================================================== |
|---|
| 3 | --- conv_temp_eden_ioniz.cpp (revision 3604) |
|---|
| 4 | +++ conv_temp_eden_ioniz.cpp (working copy) |
|---|
| 5 | @@ -96,7 +96,6 @@ |
|---|
| 6 | // here starts the standard solver for variable temperature |
|---|
| 7 | iter_track TeTrack; |
|---|
| 8 | double t1=0, error1=0, t2, error2; |
|---|
| 9 | - bool lgHighTeSearch = false; |
|---|
| 10 | |
|---|
| 11 | t2 = phycon.te; |
|---|
| 12 | error2 = CoolHeatError( t2 ); |
|---|
| 13 | @@ -109,38 +108,33 @@ |
|---|
| 14 | |
|---|
| 15 | TeTrack.clear(); |
|---|
| 16 | |
|---|
| 17 | + error1 = error2; |
|---|
| 18 | + t1 = t2; |
|---|
| 19 | + |
|---|
| 20 | + // the factor 1.2 creates 20% safety margin |
|---|
| 21 | + |
|---|
| 22 | + if ( conv.dCmHdT != 0.0) |
|---|
| 23 | + step = abs(safe_div( error1, conv.dCmHdT, 0. )); |
|---|
| 24 | + else |
|---|
| 25 | + step = 0.0; |
|---|
| 26 | + |
|---|
| 27 | // set up an initial guess for the bracket |
|---|
| 28 | // t2 was already initialized outside the main loop, or is copied from the |
|---|
| 29 | // previous iteration. don't record this evaluation, it may be poorly converged |
|---|
| 30 | - for( int i=0; i < 2 && !lgAbort; ++i ) |
|---|
| 31 | + for( int i=0; error1*error2 > 0. && i < 100 && !lgAbort; ++i ) |
|---|
| 32 | { |
|---|
| 33 | t1 = t2; |
|---|
| 34 | error1 = error2; |
|---|
| 35 | |
|---|
| 36 | - // the factor 1.2 creates 20% safety margin |
|---|
| 37 | - step = safe_div( -1.2*error1, conv.dCmHdT, 0. ); |
|---|
| 38 | - step = sign( min( abs(step), factor*t1 ), step ); |
|---|
| 39 | - t2 = t1 + step; |
|---|
| 40 | + double maxstep = factor*t1; |
|---|
| 41 | + step = 2.0*step; |
|---|
| 42 | + if (step == 0.0 || step > maxstep) |
|---|
| 43 | + step = maxstep; |
|---|
| 44 | + t2 = t1 + sign(step,-error1 ); |
|---|
| 45 | error2 = CoolHeatError( t2 ); |
|---|
| 46 | TeTrack.add( t2, error2 ); |
|---|
| 47 | } |
|---|
| 48 | |
|---|
| 49 | - int j = 0; |
|---|
| 50 | - |
|---|
| 51 | - // now hunt until we have bracketed the solution |
|---|
| 52 | - while( error1*error2 > 0. && j++ < 2*DEF_ITER && !lgAbort ) |
|---|
| 53 | - { |
|---|
| 54 | - t1 = t2; |
|---|
| 55 | - error1 = error2; |
|---|
| 56 | - double deriv = TeTrack.deriv(7); |
|---|
| 57 | - // the factor 1.2 creates 20% safety margin |
|---|
| 58 | - step = safe_div( -1.2*error1, deriv, 0. ); |
|---|
| 59 | - step = sign( min( abs(step), factor*t1 ), step ); |
|---|
| 60 | - t2 = t1 + step; |
|---|
| 61 | - error2 = CoolHeatError( t2 ); |
|---|
| 62 | - TeTrack.add( t2, error2 ); |
|---|
| 63 | - } |
|---|
| 64 | - |
|---|
| 65 | if( trace.nTrConvg >= 2 && error1*error2 > 0. && !lgAbort ) |
|---|
| 66 | { |
|---|
| 67 | fprintf( ioQQQ, " ConvTempEdenIoniz: bracket1 fails t1: %e %e t2: %e %e\n", |
|---|
| 68 | @@ -148,82 +142,8 @@ |
|---|
| 69 | TeTrack.print_history(); |
|---|
| 70 | } |
|---|
| 71 | |
|---|
| 72 | - // If we reach this point without having bracketed the solution this means |
|---|
| 73 | - // that searching using the derivate doesn't work... This typically happens |
|---|
| 74 | - // when the heating and cooling curve no longer intersect, which implies |
|---|
| 75 | - // that C-H vs Te now has a local maximum, but doesn't change sign any more. |
|---|
| 76 | - // Hence we have reached a thermal front and we need to take a discontinuous |
|---|
| 77 | - // step up or down. In the first while loop we search upwards, in the second |
|---|
| 78 | - // downwards. |
|---|
| 79 | - |
|---|
| 80 | - int nUP = 2; |
|---|
| 81 | - |
|---|
| 82 | - if( nzone > 4 && !lgHighTeSearch && !lgAbort ) |
|---|
| 83 | - { |
|---|
| 84 | - double x[3], y[3], a, b, siga, sigb; |
|---|
| 85 | - for( int i=0; i < 3; ++i ) |
|---|
| 86 | - { |
|---|
| 87 | - x[i] = double(struc.depth[nzone-4+i]); |
|---|
| 88 | - y[i] = double(struc.testr[nzone-4+i]); |
|---|
| 89 | - } |
|---|
| 90 | - linfit( 3, x, y, a, siga, b, sigb ); |
|---|
| 91 | - // if the temperature is rising significantly in the previous zones |
|---|
| 92 | - // we may have to take a big step upwards, so allow room for that... |
|---|
| 93 | - if( b > sigb ) |
|---|
| 94 | - { |
|---|
| 95 | - lgHighTeSearch = true; |
|---|
| 96 | - nUP = 25; |
|---|
| 97 | - } |
|---|
| 98 | - if( trace.nTrConvg >= 2 ) |
|---|
| 99 | - { |
|---|
| 100 | - fprintf( ioQQQ, " ConvTempEdenIoniz: temp deriv %g +/- %g\n", b, sigb ); |
|---|
| 101 | - if( lgHighTeSearch ) |
|---|
| 102 | - fprintf( ioQQQ, " ConvTempEdenIoniz: doing high Te search...\n" ); |
|---|
| 103 | - } |
|---|
| 104 | - } |
|---|
| 105 | - |
|---|
| 106 | - // hunting using the derivative failed, so first start hunting upwards |
|---|
| 107 | - // |
|---|
| 108 | - // guard against overrunning the maximum temperature limit of the code |
|---|
| 109 | - // so that we are guaranteed to get a chance to hunt downwards |
|---|
| 110 | - while( error1*error2 > 0. && j++ < (2+nUP)*DEF_ITER && !lgAbort && |
|---|
| 111 | - t2*(1.+factor) <= 0.99999*phycon.TEMP_LIMIT_HIGH ) |
|---|
| 112 | - { |
|---|
| 113 | - t1 = t2; |
|---|
| 114 | - error1 = error2; |
|---|
| 115 | - if( lgHighTeSearch && t2 > 4000. ) |
|---|
| 116 | - factor = 0.05; |
|---|
| 117 | - t2 = t1*(1.+factor); |
|---|
| 118 | - error2 = CoolHeatError( t2 ); |
|---|
| 119 | - TeTrack.add( t2, error2 ); |
|---|
| 120 | - } |
|---|
| 121 | - |
|---|
| 122 | - // hunting upwards failed, so start hunting downwards |
|---|
| 123 | - // we may need to take a big plunge downwards, so # of iter should be big |
|---|
| 124 | - // |
|---|
| 125 | - // don't guard against dropping below the minimum temperature limit of the |
|---|
| 126 | - // code so that we get an abort if the solver has completely lost the solution. |
|---|
| 127 | - // this is a serious condition that needs to be investigated... |
|---|
| 128 | - while( error1*error2 > 0. && j++ < (2+2*nUP+9)*DEF_ITER && !lgAbort ) |
|---|
| 129 | - { |
|---|
| 130 | - t1 = t2; |
|---|
| 131 | - error1 = error2; |
|---|
| 132 | - if( lgHighTeSearch && t2 <= 4000. ) |
|---|
| 133 | - factor = DEF_FACTOR; |
|---|
| 134 | - t2 = t1*(1.-factor); |
|---|
| 135 | - error2 = CoolHeatError( t2 ); |
|---|
| 136 | - TeTrack.add( t2, error2 ); |
|---|
| 137 | - } |
|---|
| 138 | - |
|---|
| 139 | factor = DEF_FACTOR; |
|---|
| 140 | |
|---|
| 141 | - if( trace.nTrConvg >= 2 && error1*error2 > 0. && !lgAbort ) |
|---|
| 142 | - { |
|---|
| 143 | - fprintf( ioQQQ, " ConvTempEdenIoniz: bracket2 fails t1: %e %e t2: %e %e\n", |
|---|
| 144 | - t1, error1, t2, error2 ); |
|---|
| 145 | - TeTrack.print_history(); |
|---|
| 146 | - } |
|---|
| 147 | - |
|---|
| 148 | // keeping the history up until now has a bad effect on convergence |
|---|
| 149 | // so we wipe the slate clean.... |
|---|
| 150 | TeTrack.clear(); |
|---|
