Changeset 2126

Show
Ignore:
Timestamp:
06/24/08 16:48:17 (2 months ago)
Author:
gary
Message:

conv_base.cpp, radius_increment.cpp move test for consistency between iso level and ionization solvers from radius_increment.cpp to conv_base.cpp, also increase error to 3% - enter ticket 70 describing this

conv_fail.cpp - improve diagnostic prints upon failure

Location:
trunk/source
Files:
3 modified

Legend:

Unmodified
Added
Removed
  • trunk/source/conv_base.cpp

    r2119 r2126  
    890890        } 
    891891 
     892        /* check that iso level population and ionization balance solvers 
     893         * are consistent - problems with this test are described in Ticket 70 */ 
     894        if( !dynamics.lgAdvection && !conv.lgSearch ) 
     895        { 
     896                for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO ) 
     897                { 
     898                        for( nelem=ipISO; nelem<LIMELM; ++nelem ) 
     899                        { 
     900                                if( dense.lgElmtOn[nelem] &&  
     901                                        dense.xIonDense[nelem][nelem-ipISO]/dense.eden>conv.EdenErrorAllowed / 10. && 
     902                                        !(iso.chTypeAtomUsed[ipISO][nelem][0]=='L' &&  
     903                                        iso.xIonSimple[ipISO][nelem]<1e-10) ) 
     904                                { 
     905                                        /* we will check that ground pops are within this factor of the xIonDense value 
     906                                        * these are only converged down to some extent - probably this should be related 
     907                                        * to the ionization convergence error */ 
     908                                        const realnum err_tol = 3e-2f; 
     909                                        realnum abund = 0.; 
     910                                        for (long n=0; n<iso.numLevels_local[ipISO][nelem]; ++n) 
     911                                        { 
     912                                                abund += (realnum)StatesElem[ipISO][nelem][n].Pop; 
     913                                        } 
     914                                        /* make sure that populations are valid, first check Pop2Ion  */ 
     915                                        if( iso.pop_ion_ov_neut[ipISO][nelem] > 1e-15 && 
     916                                                ! fp_equal_tol(abund*dense.xIonDense[nelem][nelem+1-ipISO],dense.xIonDense[nelem][nelem-ipISO], 
     917                                                err_tol*dense.xIonDense[nelem][nelem-ipISO]) 
     918                                                ) 
     919 
     920                                        if( 
     921                                                iso.pop_ion_ov_neut[ipISO][nelem] > 1e-15 && 
     922                                                fabs(abund*dense.xIonDense[nelem][nelem+1-ipISO]-dense.xIonDense[nelem][nelem-ipISO])> 
     923                                                err_tol*dense.xIonDense[nelem][nelem-ipISO] )  
     924                                        { 
     925                                                /* not converged */ 
     926                                                conv.lgConvIoniz = false; 
     927                                                strcpy( conv.chConvIoniz, "ISO solvers agree"); 
     928                                                /* indicates which species */ 
     929                                                conv.BadConvIoniz[0] = ipISO*100 + nelem; 
     930                                                conv.BadConvIoniz[1] =  
     931                                                        (abund*dense.xIonDense[nelem][nelem+1-ipISO]- 
     932                                                        dense.xIonDense[nelem][nelem-ipISO])/ 
     933                                                        SDIV(dense.xIonDense[nelem][nelem-ipISO]); 
     934                                                /*fprintf(ioQQQ,"DEBUG ISO agree %g err %.2e abund %.2e\n", 
     935                                                        conv.BadConvIoniz[0] , 
     936                                                        conv.BadConvIoniz[1], 
     937                                                        dense.xIonDense[nelem][nelem-ipISO]);*/ 
     938                                        } 
     939                                } 
     940                        } 
     941                } 
     942        } 
     943 
    892944        /* update some counters that keep track of how many times this routine 
    893945         * has been called */ 
  • trunk/source/conv_fail.cpp

    r1771 r2126  
    119119                if( called.lgTalk ) 
    120120                { 
    121                         fprintf( ioQQQ, " PROBLEM  ConvFail %li, %s ionization not converged iteration %li zone %li fnzone %.2f \n",  
     121                        fprintf( ioQQQ, " PROBLEM  ConvFail %li, %s ionization not converged" 
     122                                " iteration %li zone %li fnzone %.2f reason %s BadConvIoniz0:%g [1]=%g\n",  
    122123                          conv.nIonFail,  
    123124                          chDetail , 
    124125                          iteration , 
    125126                          nzone, 
    126                           fnzone  ); 
     127                          fnzone , 
     128                          conv.chConvIoniz, 
     129                          conv.BadConvIoniz[0], 
     130                          conv.BadConvIoniz[1]); 
    127131                } 
    128132        } 
     
    135139                if( called.lgTalk ) 
    136140                { 
    137                         fprintf( ioQQQ, " PROBLEM  ConvFail %li, %s population not converged iteration %li zone %li fnzone %.2f \n",  
     141                        fprintf( ioQQQ, " PROBLEM  ConvFail %li, %s population not converged" 
     142                                " iteration %li zone %li fnzone %.2f %s %g %g\n",  
    138143                          conv.nPopFail,  
    139144                          chDetail , 
    140145                          iteration, 
    141146                          nzone ,  
    142                           fnzone  ); 
     147                          fnzone  , 
     148                          conv.chConvIoniz, 
     149                          conv.BadConvIoniz[0], 
     150                          conv.BadConvIoniz[1]); 
    143151                } 
    144152        } 
     
    150158                if( called.lgTalk ) 
    151159                { 
    152                         fprintf( ioQQQ, " PROBLEM  ConvFail %ld, a grain failure occurred iteration %li zone %li fnzone  %.2f \n",  
     160                        fprintf( ioQQQ, " PROBLEM  ConvFail %ld, a grain failure occurred" 
     161                                " iteration %li zone %li fnzone  %.2f %s %g %g\n",  
    153162                          conv.nGrainFail,  
    154163                          iteration ,  
    155164                          nzone , 
    156                           fnzone ); 
     165                          fnzone , 
     166                          conv.chConvIoniz, 
     167                          conv.BadConvIoniz[0], 
     168                          conv.BadConvIoniz[1]); 
    157169                } 
    158170        } 
  • trunk/source/radius_increment.cpp

    r2092 r2126  
    7979        } 
    8080 
    81         /* following block only set of asserts to check that iso populations are sane */ 
    82 #       if !defined(NDEBUG) 
    83         if( !dynamics.lgAdvection ) 
    84         { 
    85                 long int ipISO; 
    86                 for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO ) 
    87                 { 
    88                         for( nelem=ipISO; nelem<LIMELM; ++nelem ) 
    89                         { 
    90                                 /* >>chng 05 feb 05, rm div by gas_phase for element, to eden,  
    91                                  * see explain for this date below */ 
    92                                 if( dense.lgElmtOn[nelem] &&  
    93                                         dense.xIonDense[nelem][nelem-ipISO]/dense.eden>conv.EdenErrorAllowed / 10. && 
    94                                         !(iso.chTypeAtomUsed[ipISO][nelem][0]=='L' &&  
    95                                         iso.xIonSimple[ipISO][nelem]<1e-10) ) 
    96                                 { 
    97                                         /* we will check that ground pops are within this factor of the xIonDense value 
    98                                          * these are only converged down to some extent - probably this should be related 
    99                                          * to the ionization convergence error */ 
    100                                         /* >>chng 05 feb 05, for very highly ionized sims, where N or O is He-like, 
    101                                          * the error can approach 1%, and depends on the convergence criteria. 
    102                                          * the code does not try to converge these quantities in an absolute sense,  
    103                                          * only the quantities relative to the electron or hydrogen density.  so it is not 
    104                                          * safe to do the assert for small values  
    105                                          * change is to only do this branch if abundance is large enough to be detected by 
    106                                          * the ionization convergence solvers */ 
    107                                         const realnum err_tol = 1e-2f; 
    108                                         /* >>chng 05 sep 02, when low-T solver used solution is approximate, 
    109                                          * and it must not matter (lot-T solver should not be used if it 
    110                                          * does matter - so use more liberal expectation */ 
    111                                         //if( iso.chTypeAtomUsed[ipISO][nelem][0]=='L' ) 
    112                                         //      err_chk *= 10.; 
    113                                         realnum abund = 0.; 
    114                                         for (long n=0; n<iso.numLevels_local[ipISO][nelem]; ++n) 
    115                                         { 
    116                                                 abund += (realnum)StatesElem[ipISO][nelem][n].Pop; 
    117                                         } 
    118                                         //fprintf(ioQQQ,"%.2f %ld %ld %11.4g %11.4g\n",fnzone,ipISO,nelem, 
    119                                         //                              abund*dense.xIonDense[nelem][nelem+1-ipISO]/SDIV(dense.xIonDense[nelem][nelem-ipISO])-1, 
    120                                         //                              dense.xIonDense[nelem][nelem-ipISO]); 
    121                                         /* make sure that populations are valid, first check Pop2Ion  */ 
    122                                         if( iso.pop_ion_ov_neut[ipISO][nelem] > 1e-15 && 
    123                                                         ! fp_equal_tol(abund*dense.xIonDense[nelem][nelem+1-ipISO],dense.xIonDense[nelem][nelem-ipISO], 
    124                                                                                                  err_tol*dense.xIonDense[nelem][nelem-ipISO]) 
    125                                                 /*>>chng 06 jul 22, only do check if pops are within tolerance of double */ 
    126                                                  ) 
    127                                         { 
    128                                                 fprintf(ioQQQ, 
    129                                                         " DISASTER radius_increment finds inconsistent populations, Pop2Ion zn %.2f ipISO %li nelem %li %.4e %.4e solver:%s\n",  
    130                                                         fnzone, 
    131                                                         ipISO , nelem , 
    132                                                         abund*dense.xIonDense[nelem][nelem+1-ipISO], 
    133                                                         dense.xIonDense[nelem][nelem-ipISO] , 
    134                                                         iso.chTypeAtomUsed[ipISO][nelem]); 
    135                                                 fprintf(ioQQQ," level solver %s found pop_ion_ov_neut of %.5e", 
    136                                                         iso.chTypeAtomUsed[ipISO][nelem] , 
    137                                                         iso.pop_ion_ov_neut[ipISO][nelem] ); 
    138                                                 fprintf(ioQQQ," ion_solver found pop_ion_ov_neut of %.5e\n", 
    139                                                         dense.xIonDense[nelem][nelem+1-ipISO]/SDIV( dense.xIonDense[nelem][nelem-ipISO] ) ); 
    140                                                 fprintf(ioQQQ, 
    141                                                         "simple %.3e Test shows Pop2Ion (%.6e) != xIonDense[nelem]/[nelem+1] (%.6e) \n",  
    142                                                         iso.xIonSimple[ipISO][nelem], 
    143                                                         abund , 
    144                                                         dense.xIonDense[nelem][nelem-ipISO]/(double)SDIV(dense.xIonDense[nelem][nelem+1-ipISO]) ); 
    145                                         } 
    146  
    147                                         /* >>chng 06 jul 24, add assert that this is not search  
    148                                          * phase to confirm that removing other asserts on this  
    149                                          * are ok - we cannot be in search phase in this routine  
    150                                          * if ok, rm this and all ref to lgSearch, perhaps conv.h header */ 
    151                                         ASSERT( !conv.lgSearch ); 
    152                                         ASSERT( /*>>chng 06 jul 22, only do check if pops are within tolerance of double */ 
    153                                                 iso.pop_ion_ov_neut[ipISO][nelem] < 1e-15 || 
    154                                                 fp_equal_tol(abund*dense.xIonDense[nelem][nelem+1-ipISO],dense.xIonDense[nelem][nelem-ipISO], 
    155                                                                                                  err_tol*dense.xIonDense[nelem][nelem-ipISO]) ); 
    156                                 } 
    157                         } 
    158                 } 
    159         } 
    160 #       endif 
    161  
    16281        if( trace.lgTrace ) 
    16382        {